diff --git a/libs/solver/direct.h b/libs/solver/direct.h index a5fb8ef..7a2aa22 100644 --- a/libs/solver/direct.h +++ b/libs/solver/direct.h @@ -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& a_matrix); /** * @brief Solves the linear system, using the solver's internally stored factorised coefficient matrix. @@ -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& x_vector, const Vector_Dense& b_vector){ return static_cast<_solver*>(this)->solve_system(x_vector, b_vector); }; diff --git a/libs/solver/direct_lower_upper_factorisation.h b/libs/solver/direct_lower_upper_factorisation.h index bce4969..ea4d796 100644 --- a/libs/solver/direct_lower_upper_factorisation.h +++ b/libs/solver/direct_lower_upper_factorisation.h @@ -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& a_matrix); /** * @brief Solves the linear system, using the solver's internally stored factorised coefficient matrix. @@ -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& x_vector, const Vector_Dense& b_vector); /** * @brief Gets the current configuration of the solver. @@ -104,7 +104,7 @@ class Direct_Lower_Upper_Factorisation private: bool factorised{false}; // lu_factorised; // lu_factorised; // pivots; // bool Direct_Lower_Upper_Factorisation<_solver_type, _size, _pivot>::factorise( - const Matrix_Dense<_size, _size>& a_matrix) { + const Matrix_Dense& a_matrix) { // Initialise factorisation data. factorised = false; @@ -105,7 +105,7 @@ bool Direct_Lower_Upper_Factorisation<_solver_type, _size, _pivot>::factorise( */ template 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& x_vector, const Vector_Dense& b_vector) { ASSERT_DEBUG(b_vector.size() == lu_factorised.size_row(), "Constant vector not of the correct size."); diff --git a/libs/solver/solver.h b/libs/solver/solver.h index 6d61ff3..7fd1bf9 100644 --- a/libs/solver/solver.h +++ b/libs/solver/solver.h @@ -33,13 +33,13 @@ namespace Disa{ // Forward declarations class Matrix_Sparse; -template class Vector_Dense; +template 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: @@ -53,8 +53,8 @@ class Solver { std::unique_ptr, 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& x_vector, + const Vector_Dense& b_vector) { switch(solver.index()) { case 0: ERROR("LU solver does not support sparse matrices."); @@ -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& x_vector, const Vector_Dense& b_vector) { switch(solver.index()) { case 0: return std::get > >(solver)->solve_system(x_vector, b_vector); diff --git a/libs/solver/solver_fixed_point.cpp b/libs/solver/solver_fixed_point.cpp index 247cc33..16bbda0 100644 --- a/libs/solver/solver_fixed_point.cpp +++ b/libs/solver/solver_fixed_point.cpp @@ -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& x_vector, Vector_Dense& x_update, + const Vector_Dense& b_vector, const Scalar omega = 1) { // forward sweep FOR(i_row, a_matrix.size_row()) { Scalar offs_row_dot = 0; @@ -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& x_vector, Vector_Dense& x_update, + const Vector_Dense& b_vector, const Scalar omega = 1){ // forward sweep for(auto i_row = a_matrix.size_row() - 1; i_row != std::numeric_limits::max(); --i_row) { Scalar offs_row_dot = 0; @@ -61,9 +61,10 @@ void Solver_Fixed_Point::in } template<> -Convergence_Data Solver_Fixed_Point::solve_system(const Matrix_Sparse& a_matrix, - Vector_Dense<0>& x_vector, - const Vector_Dense<0>& b_vector) { +Convergence_Data +Solver_Fixed_Point::solve_system(const Matrix_Sparse& a_matrix, + Vector_Dense& x_vector, + const Vector_Dense& b_vector) { data.working.resize(a_matrix.size_row()); Convergence_Data convergence_data = Convergence_Data(); while(!data.limits.is_converged(convergence_data)) { @@ -83,9 +84,10 @@ void Solver_Fixed_Point::ini } template<> -Convergence_Data Solver_Fixed_Point::solve_system(const Matrix_Sparse& a_matrix, - Vector_Dense<0>& x_vector, - const Vector_Dense<0>& b_vector) { +Convergence_Data +Solver_Fixed_Point::solve_system(const Matrix_Sparse& a_matrix, + Vector_Dense& x_vector, + const Vector_Dense& 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); @@ -106,7 +108,7 @@ void Solver_Fixed_Point Convergence_Data Solver_Fixed_Point::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& x_vector, const Vector_Dense& 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); diff --git a/libs/solver/solver_fixed_point.h b/libs/solver/solver_fixed_point.h index 0af9e71..94f1fae 100644 --- a/libs/solver/solver_fixed_point.h +++ b/libs/solver/solver_fixed_point.h @@ -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 working; }; struct Solver_Fixed_Point_Sor_Data : public Solver_Data{ @@ -58,8 +58,8 @@ class Solver_Fixed_Point : public Solver_Iterative& x_vector, - const Vector_Dense<0>& b_vector); + Convergence_Data solve_system(const Matrix_Sparse& a_matrix, Vector_Dense& x_vector, + const Vector_Dense& b_vector); }; typedef Solver_Fixed_Point Solver_Jacobi; diff --git a/libs/solver/solver_iterative.h b/libs/solver/solver_iterative.h index 89789b6..82572d7 100644 --- a/libs/solver/solver_iterative.h +++ b/libs/solver/solver_iterative.h @@ -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& x_vector, + const Vector_Dense& b_vector){ return static_cast<_solver*>(this)->solve_system(matrix, x_vector, b_vector); }; diff --git a/libs/solver/tests/test_direct.cpp b/libs/solver/tests/test_direct.cpp index 1792d79..f24e55a 100644 --- a/libs/solver/tests/test_direct.cpp +++ b/libs/solver/tests/test_direct.cpp @@ -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 b_vector = {6, 2}; + Vector_Dense x_vector = {0, 0, 0}; + Matrix_Dense 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), "./*"); @@ -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 b_vector = {6, 2, 7}; + Vector_Dense x_vector = {0, 0, 0}; + Matrix_Dense matrix = {{0, 0, 0}, {9, 5, 1}, {4, 3, 8}}; // singular. EXPECT_FALSE(lu_solver.factorise(matrix)); EXPECT_FALSE(lup_solver.factorise(matrix)); @@ -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 matrix = {{2, 7, 6}, {9, 5, 1}, {4, 3, 8}}; + Vector_Dense b_vector = {6, 2, 7}; + Vector_Dense x_vector = {0, 0, 0}; + Vector_Dense solution = {0.04166666666666696, 0.16666666666666657, 0.7916666666666667}; // Test LU EXPECT_TRUE(lu_solver.factorise(matrix)); diff --git a/libs/solver/tests/test_solver.cpp b/libs/solver/tests/test_solver.cpp index bd8cff7..c8a176f 100644 --- a/libs/solver/tests/test_solver.cpp +++ b/libs/solver/tests/test_solver.cpp @@ -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 x_vector; //!< Solution vector of the linear system. + Vector_Dense b_vector; //!< Constant vector of the linear system. + Vector_Dense b_vector_0; //!< Constant vector of the linear system. + Matrix_Dense 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). diff --git a/libs/solver/tests/test_solver_utilities.cpp b/libs/solver/tests/test_solver_utilities.cpp index dde9993..8bd30d1 100644 --- a/libs/solver/tests/test_solver_utilities.cpp +++ b/libs/solver/tests/test_solver_utilities.cpp @@ -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 solution; solution.resize(5, 2.0); - Vector_Dense<0> constant = {3.0, -1.0, 3.0, -1.0, 3.0}; + Vector_Dense 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); @@ -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({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); @@ -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 solution; solution.resize(5, 2.0); - Vector_Dense<0> constant = {3.0, -1.0, 3.0, -1.0, 3.0}; + Vector_Dense 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({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{}, constant), ".*"); + EXPECT_DEATH(compute_residual(Matrix_Sparse(1, 1), solution, Vector_Dense{}), ".*"); } \ No newline at end of file