Skip to content

Commit

Permalink
solver:
Browse files Browse the repository at this point in the history 
- Vector and Matrix dense now use Scalar as the template type implementation.
  • Loading branch information
bevan committed Mar 9, 2024
1 parent 48613ab commit 4f4c7f5
Show file tree
Hide file tree
Showing 10 changed files with 52 additions and 53 deletions.
4 changes: 2 additions & 2 deletions libs/solver/direct.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ class Direct
* @param[in] a_matrix The coefficient matrix to factorise.
* @return True if the matrix was factorised successfully, else false.
*/
bool factorise(const Matrix_Dense<_size, _size>& a_matrix);
bool factorise(const Matrix_Dense<Scalar, _size, _size>& a_matrix);

/**
* @brief Solves the linear system, using the solver's internally stored factorised coefficient matrix.
Expand All @@ -85,7 +85,7 @@ class Direct
* @warning The factorise function must have been called before this function, and returned true for a correct solver
* to occur.
*/
const Convergence_Data& solve(Vector_Dense<_size>& x_vector, const Vector_Dense<_size>& b_vector){
const Convergence_Data& solve(Vector_Dense<Scalar, _size>& x_vector, const Vector_Dense<Scalar, _size>& b_vector){
return static_cast<_solver*>(this)->solve_system(x_vector, b_vector);
};

Expand Down
6 changes: 3 additions & 3 deletions libs/solver/direct_lower_upper_factorisation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ class Direct_Lower_Upper_Factorisation
* @param[in] a_matrix The coefficient matrix to factorise.
* @return True if the matrix was factorised successfully, else false for degenerate/singular matrices.
*/
bool factorise(const Matrix_Dense<_size, _size>& a_matrix);
bool factorise(const Matrix_Dense<Scalar, _size, _size>& a_matrix);

/**
* @brief Solves the linear system, using the solver's internally stored factorised coefficient matrix.
Expand All @@ -87,7 +87,7 @@ class Direct_Lower_Upper_Factorisation
* @warning The factorise function must have been called before this function, and returned true for a correct solver
* to occur.
*/
Convergence_Data solve_system(Vector_Dense<_size>& x_vector, const Vector_Dense<_size>& b_vector);
Convergence_Data solve_system(Vector_Dense<Scalar, _size>& x_vector, const Vector_Dense<Scalar, _size>& b_vector);

/**
* @brief Gets the current configuration of the solver.
Expand All @@ -104,7 +104,7 @@ class Direct_Lower_Upper_Factorisation
private:
bool factorised{false}; //<! Has factorisation been completed successfully.
Scalar factorisation_tolerance{default_absolute}; //<! The value below which diagonal entires should be considered zero.
Matrix_Dense<_size, _size> lu_factorised; //<! The factorised LU coefficient matrix (implicit 1's on diagonal of L).
Matrix_Dense<Scalar, _size, _size> lu_factorised; //<! The factorised LU coefficient matrix (implicit 1's on diagonal of L).
std::vector<std::size_t> pivots; //<! The pivots indicies, is only sized for LUP solver.
};

Expand Down
4 changes: 2 additions & 2 deletions libs/solver/direct_lower_upper_factorisation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ namespace Disa {
*/
template<Solver_Type _solver_type, std::size_t _size, bool _pivot>
bool Direct_Lower_Upper_Factorisation<_solver_type, _size, _pivot>::factorise(
const Matrix_Dense<_size, _size>& a_matrix) {
const Matrix_Dense<Scalar, _size, _size>& a_matrix) {

// Initialise factorisation data.
factorised = false;
Expand Down Expand Up @@ -105,7 +105,7 @@ bool Direct_Lower_Upper_Factorisation<_solver_type, _size, _pivot>::factorise(
*/
template<Solver_Type _solver_type, std::size_t _size, bool _pivot>
Convergence_Data Direct_Lower_Upper_Factorisation<_solver_type, _size, _pivot>::solve_system(
Vector_Dense<_size>& x_vector, const Vector_Dense<_size>& b_vector) {
Vector_Dense<Scalar, _size>& x_vector, const Vector_Dense<Scalar, _size>& b_vector) {

ASSERT_DEBUG(b_vector.size() == lu_factorised.size_row(), "Constant vector not of the correct size.");

Expand Down
10 changes: 5 additions & 5 deletions libs/solver/solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,13 +33,13 @@ namespace Disa{

// Forward declarations
class Matrix_Sparse;
template<std::size_t> class Vector_Dense;
template<typename, std::size_t> class Vector_Dense;

/**
* @brief Solver
* @brief General Solver Class for all solver types.
*
* This class serves at a 'all in one' solver class for both sparse and dense systems. To acheive
* This class serves at a 'all in one' solver class for both sparse and dense systems. To achieve
*/
class Solver {
public:
Expand All @@ -53,8 +53,8 @@ class Solver {
std::unique_ptr<Sover_Sor>,
std::nullptr_t> solver{nullptr};

Convergence_Data solve(const Matrix_Sparse& a_matrix, Vector_Dense<0>& x_vector,
const Vector_Dense<0>& b_vector) {
Convergence_Data solve(const Matrix_Sparse& a_matrix, Vector_Dense<Scalar, 0>& x_vector,
const Vector_Dense<Scalar, 0>& b_vector) {

switch(solver.index()) {
case 0: ERROR("LU solver does not support sparse matrices.");
Expand All @@ -68,7 +68,7 @@ class Solver {
}
};

Convergence_Data solve(Vector_Dense<0>& x_vector, const Vector_Dense<0>& b_vector) {
Convergence_Data solve(Vector_Dense<Scalar, 0>& x_vector, const Vector_Dense<Scalar, 0>& b_vector) {

switch(solver.index()) {
case 0: return std::get<std::unique_ptr<Solver_LU<0> > >(solver)->solve_system(x_vector, b_vector);
Expand Down
24 changes: 13 additions & 11 deletions libs/solver/solver_fixed_point.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@
namespace Disa {

inline void forward_sweep(const Matrix_Sparse& a_matrix,
const Vector_Dense<0>& x_vector, Vector_Dense<0>& x_update,
const Vector_Dense<0>& b_vector, const Scalar omega = 1) {
const Vector_Dense<Scalar, 0>& x_vector, Vector_Dense<Scalar, 0>& x_update,
const Vector_Dense<Scalar, 0>& b_vector, const Scalar omega = 1) {
// forward sweep
FOR(i_row, a_matrix.size_row()) {
Scalar offs_row_dot = 0;
Expand All @@ -40,8 +40,8 @@ inline void forward_sweep(const Matrix_Sparse& a_matrix,
}

inline void backward_sweep(const Matrix_Sparse& a_matrix,
const Vector_Dense<0>& x_vector, Vector_Dense<0>& x_update,
const Vector_Dense<0>& b_vector, const Scalar omega = 1){
const Vector_Dense<Scalar, 0>& x_vector, Vector_Dense<Scalar, 0>& x_update,
const Vector_Dense<Scalar, 0>& b_vector, const Scalar omega = 1){
// forward sweep
for(auto i_row = a_matrix.size_row() - 1; i_row != std::numeric_limits<std::size_t>::max(); --i_row) {
Scalar offs_row_dot = 0;
Expand All @@ -61,9 +61,10 @@ void Solver_Fixed_Point<Solver_Type::jacobi, Solver_Fixed_Point_Jacobi_Data>::in
}

template<>
Convergence_Data Solver_Fixed_Point<Solver_Type::jacobi, Solver_Fixed_Point_Jacobi_Data>::solve_system(const Matrix_Sparse& a_matrix,
Vector_Dense<0>& x_vector,
const Vector_Dense<0>& b_vector) {
Convergence_Data
Solver_Fixed_Point<Solver_Type::jacobi, Solver_Fixed_Point_Jacobi_Data>::solve_system(const Matrix_Sparse& a_matrix,
Vector_Dense<Scalar, 0>& x_vector,
const Vector_Dense<Scalar, 0>& b_vector) {
data.working.resize(a_matrix.size_row());
Convergence_Data convergence_data = Convergence_Data();
while(!data.limits.is_converged(convergence_data)) {
Expand All @@ -83,9 +84,10 @@ void Solver_Fixed_Point<Solver_Type::gauss_seidel, Solver_Fixed_Point_Data>::ini
}

template<>
Convergence_Data Solver_Fixed_Point<Solver_Type::gauss_seidel, Solver_Fixed_Point_Data>::solve_system(const Matrix_Sparse& a_matrix,
Vector_Dense<0>& x_vector,
const Vector_Dense<0>& b_vector) {
Convergence_Data
Solver_Fixed_Point<Solver_Type::gauss_seidel, Solver_Fixed_Point_Data>::solve_system(const Matrix_Sparse& a_matrix,
Vector_Dense<Scalar, 0>& x_vector,
const Vector_Dense<Scalar, 0>& b_vector) {
Convergence_Data convergence_data = Convergence_Data();
while(!data.limits.is_converged(convergence_data)) {
forward_sweep(a_matrix, x_vector, x_vector, b_vector, 1.0);
Expand All @@ -106,7 +108,7 @@ void Solver_Fixed_Point<Solver_Type::successive_over_relaxation, Solver_Fixed_Po

template<>
Convergence_Data Solver_Fixed_Point<Solver_Type::successive_over_relaxation, Solver_Fixed_Point_Sor_Data>::solve_system(
const Matrix_Sparse& a_matrix, Vector_Dense<0>& x_vector, const Vector_Dense<0>& b_vector) {
const Matrix_Sparse& a_matrix, Vector_Dense<Scalar, 0>& x_vector, const Vector_Dense<Scalar, 0>& b_vector) {
Convergence_Data convergence_data = Convergence_Data();
while(!data.limits.is_converged(convergence_data)) {
forward_sweep(a_matrix, x_vector, x_vector, b_vector, 1.5);
Expand Down
6 changes: 3 additions & 3 deletions libs/solver/solver_fixed_point.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ struct Solver_Fixed_Point_Data : public Solver_Data{
};

struct Solver_Fixed_Point_Jacobi_Data : public Solver_Data{
Vector_Dense<0> working;
Vector_Dense<Scalar, 0> working;
};

struct Solver_Fixed_Point_Sor_Data : public Solver_Data{
Expand Down Expand Up @@ -58,8 +58,8 @@ class Solver_Fixed_Point : public Solver_Iterative<Solver_Fixed_Point<_solver_ty
* @param b_vector
* @return
*/
Convergence_Data solve_system(const Matrix_Sparse& a_matrix, Vector_Dense<0>& x_vector,
const Vector_Dense<0>& b_vector);
Convergence_Data solve_system(const Matrix_Sparse& a_matrix, Vector_Dense<Scalar, 0>& x_vector,
const Vector_Dense<Scalar, 0>& b_vector);
};

typedef Solver_Fixed_Point<Solver_Type::jacobi, Solver_Fixed_Point_Jacobi_Data> Solver_Jacobi;
Expand Down
4 changes: 2 additions & 2 deletions libs/solver/solver_iterative.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,8 @@ class Solver_Iterative
return static_cast<_solver*>(this)->initialise_solver(solver_config);
};

const Convergence_Data& solve(const Matrix_Sparse& matrix, Vector_Dense<0>& x_vector,
const Vector_Dense<0>& b_vector){
const Convergence_Data& solve(const Matrix_Sparse& matrix, Vector_Dense<Scalar, 0>& x_vector,
const Vector_Dense<Scalar, 0>& b_vector){
return static_cast<_solver*>(this)->solve_system(matrix, x_vector, b_vector);
};

Expand Down
20 changes: 10 additions & 10 deletions libs/solver/tests/test_direct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,9 @@ TEST_F(direct_solvers, lower_upper_factorisation_death_test) {
EXPECT_DEATH(Solver_LUP<0>lup(data), "./*");

// incorrect sizes
Vector_Dense<0> b_vector = {6, 2};
Vector_Dense<0> x_vector = {0, 0, 0};
Matrix_Dense<0, 0> matrix = {{2, 7, 6}, {9, 5, 1}, {4, 3, 8}};
Vector_Dense<Scalar, 0> b_vector = {6, 2};
Vector_Dense<Scalar, 0> x_vector = {0, 0, 0};
Matrix_Dense<Scalar, 0, 0> matrix = {{2, 7, 6}, {9, 5, 1}, {4, 3, 8}};
lu_solver.factorise(matrix);
lup_solver.factorise(matrix);
EXPECT_DEATH(lu_solver.solve_system(x_vector, b_vector), "./*");
Expand All @@ -96,9 +96,9 @@ TEST_F(direct_solvers, lower_upper_factorisation_death_test) {

TEST_F(direct_solvers, lower_upper_factorisation_not_factorised) {

Vector_Dense<0> b_vector = {6, 2, 7};
Vector_Dense<0> x_vector = {0, 0, 0};
Matrix_Dense<0, 0> matrix = {{0, 0, 0}, {9, 5, 1}, {4, 3, 8}}; // singular.
Vector_Dense<Scalar, 0> b_vector = {6, 2, 7};
Vector_Dense<Scalar, 0> x_vector = {0, 0, 0};
Matrix_Dense<Scalar, 0, 0> matrix = {{0, 0, 0}, {9, 5, 1}, {4, 3, 8}}; // singular.

EXPECT_FALSE(lu_solver.factorise(matrix));
EXPECT_FALSE(lup_solver.factorise(matrix));
Expand All @@ -110,10 +110,10 @@ TEST_F(direct_solvers, lower_upper_factorisation_not_factorised) {
TEST_F(direct_solvers, lower_upper_factorise_solve) {

// Setup linear system (source https://lampx.tugraz.at/~hadley/num/ch2/2.3a.php)
Matrix_Dense<0, 0> matrix = {{2, 7, 6}, {9, 5, 1}, {4, 3, 8}};
Vector_Dense<0> b_vector = {6, 2, 7};
Vector_Dense<0> x_vector = {0, 0, 0};
Vector_Dense<0> solution = {0.04166666666666696, 0.16666666666666657, 0.7916666666666667};
Matrix_Dense<Scalar, 0, 0> matrix = {{2, 7, 6}, {9, 5, 1}, {4, 3, 8}};
Vector_Dense<Scalar, 0> b_vector = {6, 2, 7};
Vector_Dense<Scalar, 0> x_vector = {0, 0, 0};
Vector_Dense<Scalar, 0> solution = {0.04166666666666696, 0.16666666666666657, 0.7916666666666667};

// Test LU
EXPECT_TRUE(lu_solver.factorise(matrix));
Expand Down
11 changes: 4 additions & 7 deletions libs/solver/tests/test_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,10 @@ class Laplace2DProblem : public ::testing::Test
public:
Matrix_Sparse a_sparse; //!< Sparse coefficient matrix of the linear system
Matrix_Sparse a_sparse_0; //!< Sparse coefficient matrix of the linear system
Vector_Dense<0> x_vector; //!< Solution vector of the linear system.
Vector_Dense<0> b_vector; //!< Constant vector of the linear system.
Vector_Dense<0> b_vector_0; //!< Constant vector of the linear system.



Matrix_Dense<0,0> a_dense; //!< Dense coefficient matrix of the linear system
Vector_Dense<Scalar, 0> x_vector; //!< Solution vector of the linear system.
Vector_Dense<Scalar, 0> b_vector; //!< Constant vector of the linear system.
Vector_Dense<Scalar, 0> b_vector_0; //!< Constant vector of the linear system.
Matrix_Dense<Scalar, 0,0> a_dense; //!< Dense coefficient matrix of the linear system

/**
* @brief Calls below setup function with a mesh size of 100 grid points (10 on each axis).
Expand Down
16 changes: 8 additions & 8 deletions libs/solver/tests/test_solver_utilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,9 +40,9 @@ TEST(solver_utilites, update_convergence) {
matrix_sparse[2][2] = 0.5;
matrix_sparse[3][3] = 0.5;
matrix_sparse[4][4] = 0.5;
Vector_Dense<0> solution;
Vector_Dense<Scalar, 0> solution;
solution.resize(5, 2.0);
Vector_Dense<0> constant = {3.0, -1.0, 3.0, -1.0, 3.0};
Vector_Dense<Scalar, 0> constant = {3.0, -1.0, 3.0, -1.0, 3.0};
auto [weighted_l2_norm_0, l_inf_norm_0] = compute_residual(matrix_sparse, solution, constant);


Expand All @@ -59,7 +59,7 @@ TEST(solver_utilites, update_convergence) {
EXPECT_EQ(data.residual_max_normalised, 1.0);

// Second iteration, check that which must be correct after instantiation.
constant = Vector_Dense<0>({2.0, 0.0, 3.0, 0.0, 2.0});
constant = Vector_Dense<Scalar, 0>({2.0, 0.0, 3.0, 0.0, 2.0});
auto [weighted_l2_norm, l_inf_norm] = compute_residual(matrix_sparse, solution, constant);

data.update(matrix_sparse, solution, constant);
Expand Down Expand Up @@ -111,21 +111,21 @@ TEST(solver_utilites, compute_residual) {
matrix_sparse[2][2] = 0.5;
matrix_sparse[3][3] = 0.5;
matrix_sparse[4][4] = 0.5;
Vector_Dense<0> solution;
Vector_Dense<Scalar, 0> solution;
solution.resize(5, 2.0);

Vector_Dense<0> constant = {3.0, -1.0, 3.0, -1.0, 3.0};
Vector_Dense<Scalar, 0> constant = {3.0, -1.0, 3.0, -1.0, 3.0};
auto [weighted_l2_norm, l_inf_norm] = compute_residual(matrix_sparse, solution, constant);
EXPECT_DOUBLE_EQ(weighted_l2_norm, 2.0);
EXPECT_DOUBLE_EQ(l_inf_norm, 2.0);

constant = Vector_Dense<0>({2.0, 0.0, 3.0, 0.0, 2.0});
constant = Vector_Dense<Scalar, 0>({2.0, 0.0, 3.0, 0.0, 2.0});
std::tie(weighted_l2_norm, l_inf_norm) = compute_residual(matrix_sparse, solution, constant);
EXPECT_DOUBLE_EQ(weighted_l2_norm, std::sqrt(8.0/5.0));
EXPECT_DOUBLE_EQ(l_inf_norm, 2.0);

EXPECT_DEATH(compute_residual(Matrix_Sparse(5, 2), solution, constant), ".*");
EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 5), solution, constant), ".*");
EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 1), Vector_Dense<0>{}, constant), ".*");
EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 1), solution, Vector_Dense<0>{}), ".*");
EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 1), Vector_Dense<Scalar, 0>{}, constant), ".*");
EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 1), solution, Vector_Dense<Scalar, 0>{}), ".*");
}

0 comments on commit 4f4c7f5

Please sign in to comment.