diff --git a/include/LinearAlgebra/diagonalSolver.h b/include/LinearAlgebra/diagonalSolver.h index a0a92fef..f003397b 100644 --- a/include/LinearAlgebra/diagonalSolver.h +++ b/include/LinearAlgebra/diagonalSolver.h @@ -22,14 +22,8 @@ class DiagonalSolver { public: DiagonalSolver(); - DiagonalSolver(const DiagonalSolver& other); - DiagonalSolver(DiagonalSolver&& other) noexcept; - explicit DiagonalSolver(const int matrix_dimension); - DiagonalSolver& operator=(const DiagonalSolver& other); - DiagonalSolver& operator=(DiagonalSolver&& other) noexcept; - int rows() const; int columns() const; @@ -62,99 +56,45 @@ std::ostream& operator<<(std::ostream& stream, const DiagonalSolver& solver) return stream; } -// default construction template DiagonalSolver::DiagonalSolver() : matrix_dimension_(0) { } -// copy construction -template -DiagonalSolver::DiagonalSolver(const DiagonalSolver& other) - : matrix_dimension_(other.matrix_dimension_) - , diagonal_values_("Diagonal matrix values", matrix_dimension_) -{ - Kokkos::deep_copy(diagonal_values_, other.diagonal_values_); -} - -// copy assignment -template -DiagonalSolver& DiagonalSolver::operator=(const DiagonalSolver& other) -{ - if (this == &other) { - // Self-assignment, no work needed - return *this; - } - // Only allocate new memory if the sizes are different - if (matrix_dimension_ != other.matrix_dimension_) { - matrix_dimension_ = other.matrix_dimension_; - diagonal_values_ = Vector("Diagonal matrix values", matrix_dimension_); - } - Kokkos::deep_copy(diagonal_values_, other.diagonal_values_); - return *this; -} - -// move construction -template -DiagonalSolver::DiagonalSolver(DiagonalSolver&& other) noexcept - : matrix_dimension_(other.matrix_dimension_) - , diagonal_values_(std::move(other.diagonal_values_)) -{ - other.matrix_dimension_ = 0; -} - -// move assignment -template -DiagonalSolver& DiagonalSolver::operator=(DiagonalSolver&& other) noexcept -{ - matrix_dimension_ = other.matrix_dimension_; - diagonal_values_ = std::move(other.diagonal_values_); - other.matrix_dimension_ = 0; - return *this; -} - template DiagonalSolver::DiagonalSolver(const int matrix_dimension) : matrix_dimension_(matrix_dimension) , diagonal_values_("Diagonal matrix values", matrix_dimension_) { - assert(matrix_dimension_ >= 1); + assert(matrix_dimension_ >= 0); Kokkos::deep_copy(diagonal_values_, T(0)); } template int DiagonalSolver::rows() const { - assert(matrix_dimension_ >= 0); return matrix_dimension_; } template int DiagonalSolver::columns() const { - assert(matrix_dimension_ >= 0); return matrix_dimension_; } template const T& DiagonalSolver::diagonal(const int index) const { - assert(index >= 0); - assert(index < matrix_dimension_); + assert(0 <= index && index < matrix_dimension_); return diagonal_values_(index); } template T& DiagonalSolver::diagonal(const int index) { - assert(index >= 0); - assert(index < matrix_dimension_); + assert(0 <= index && index < matrix_dimension_); return diagonal_values_(index); } -// --------------- // -// Diagonal Solver // -// --------------- // - template void DiagonalSolver::solveInPlace(T* sol_rhs) const { diff --git a/include/LinearAlgebra/symmetricTridiagonalSolver.h b/include/LinearAlgebra/symmetricTridiagonalSolver.h index 64d0935e..3c0c9016 100644 --- a/include/LinearAlgebra/symmetricTridiagonalSolver.h +++ b/include/LinearAlgebra/symmetricTridiagonalSolver.h @@ -48,14 +48,8 @@ class SymmetricTridiagonalSolver { public: SymmetricTridiagonalSolver(); - SymmetricTridiagonalSolver(const SymmetricTridiagonalSolver& other); - SymmetricTridiagonalSolver(SymmetricTridiagonalSolver&& other) noexcept; - explicit SymmetricTridiagonalSolver(const int matrix_dimension); - SymmetricTridiagonalSolver& operator=(const SymmetricTridiagonalSolver& other); - SymmetricTridiagonalSolver& operator=(SymmetricTridiagonalSolver&& other) noexcept; - void is_cyclic(bool value); bool is_cyclic() const; @@ -81,11 +75,10 @@ class SymmetricTridiagonalSolver int matrix_dimension_; std::unique_ptr main_diagonal_values_; std::unique_ptr sub_diagonal_values_; - T cyclic_corner_element_ = 0.0; - bool is_cyclic_ = true; - - bool factorized_ = false; - T gamma_ = 0.0; // Used in Shermann-Morrison factorization A = B + u*v^T + T cyclic_corner_element_; + bool is_cyclic_; + bool factorized_; + T gamma_; // Used in Shermann-Morrison factorization A = B + u*v^T // Solve methods: // The notation 'u' and 'scratch' is directly taken from the implementation @@ -156,76 +149,11 @@ SymmetricTridiagonalSolver::SymmetricTridiagonalSolver() , sub_diagonal_values_(nullptr) , cyclic_corner_element_(0.0) , is_cyclic_(true) + , factorized_(false) + , gamma_(0.0) { } -// copy construction -template -SymmetricTridiagonalSolver::SymmetricTridiagonalSolver(const SymmetricTridiagonalSolver& other) - : matrix_dimension_(other.matrix_dimension_) - , main_diagonal_values_(std::make_unique(matrix_dimension_)) - , sub_diagonal_values_(std::make_unique(matrix_dimension_ - 1)) - , cyclic_corner_element_(other.cyclic_corner_element_) - , is_cyclic_(other.is_cyclic_) -{ - std::copy(other.main_diagonal_values_.get(), other.main_diagonal_values_.get() + matrix_dimension_, - main_diagonal_values_.get()); - std::copy(other.sub_diagonal_values_.get(), other.sub_diagonal_values_.get() + matrix_dimension_ - 1, - sub_diagonal_values_.get()); -} - -// copy assignment -template -SymmetricTridiagonalSolver& SymmetricTridiagonalSolver::operator=(const SymmetricTridiagonalSolver& other) -{ - if (this == &other) { - // Self-assignment, no work needed - return *this; - } - // Only allocate new memory if the sizes are different - if (matrix_dimension_ != other.matrix_dimension_) { - matrix_dimension_ = other.matrix_dimension_; - main_diagonal_values_ = std::make_unique(matrix_dimension_); - sub_diagonal_values_ = std::make_unique(matrix_dimension_ - 1); - } - cyclic_corner_element_ = other.cyclic_corner_element_; - is_cyclic_ = other.is_cyclic_; - std::copy(other.main_diagonal_values_.get(), other.main_diagonal_values_.get() + matrix_dimension_, - main_diagonal_values_.get()); - std::copy(other.sub_diagonal_values_.get(), other.sub_diagonal_values_.get() + matrix_dimension_ - 1, - sub_diagonal_values_.get()); - return *this; -} - -// move construction -template -SymmetricTridiagonalSolver::SymmetricTridiagonalSolver(SymmetricTridiagonalSolver&& other) noexcept - : matrix_dimension_(other.matrix_dimension_) - , main_diagonal_values_(std::move(other.main_diagonal_values_)) - , sub_diagonal_values_(std::move(other.sub_diagonal_values_)) - , cyclic_corner_element_(other.cyclic_corner_element_) - , is_cyclic_(other.is_cyclic_) -{ - other.matrix_dimension_ = 0; - other.cyclic_corner_element_ = 0.0; - other.is_cyclic_ = true; -} - -// move assignment -template -SymmetricTridiagonalSolver& SymmetricTridiagonalSolver::operator=(SymmetricTridiagonalSolver&& other) noexcept -{ - matrix_dimension_ = other.matrix_dimension_; - main_diagonal_values_ = std::move(other.main_diagonal_values_); - sub_diagonal_values_ = std::move(other.sub_diagonal_values_); - cyclic_corner_element_ = other.cyclic_corner_element_; - is_cyclic_ = other.is_cyclic_; - other.matrix_dimension_ = 0; - other.cyclic_corner_element_ = 0.0; - other.is_cyclic_ = true; - return *this; -} - template SymmetricTridiagonalSolver::SymmetricTridiagonalSolver(const int matrix_dimension) : matrix_dimension_(matrix_dimension) @@ -233,6 +161,8 @@ SymmetricTridiagonalSolver::SymmetricTridiagonalSolver(const int matrix_dimen , sub_diagonal_values_(std::make_unique(matrix_dimension_ - 1)) , cyclic_corner_element_(0.0) , is_cyclic_(true) + , factorized_(false) + , gamma_(0.0) { assert(matrix_dimension_ >= 1); std::fill(main_diagonal_values_.get(), main_diagonal_values_.get() + matrix_dimension_, T(0)); @@ -253,57 +183,51 @@ bool SymmetricTridiagonalSolver::is_cyclic() const template int SymmetricTridiagonalSolver::rows() const { - assert(this->matrix_dimension_ >= 0); - return this->matrix_dimension_; + return matrix_dimension_; } template int SymmetricTridiagonalSolver::columns() const { - assert(this->matrix_dimension_ >= 0); - return this->matrix_dimension_; + return matrix_dimension_; } template const T& SymmetricTridiagonalSolver::main_diagonal(const int index) const { - assert(index >= 0); - assert(index < this->matrix_dimension_); - return this->main_diagonal_values_[index]; + assert(0 <= index && index < matrix_dimension_); + return main_diagonal_values_[index]; } template T& SymmetricTridiagonalSolver::main_diagonal(const int index) { - assert(index >= 0); - assert(index < this->matrix_dimension_); - return this->main_diagonal_values_[index]; + assert(0 <= index && index < matrix_dimension_); + return main_diagonal_values_[index]; } template const T& SymmetricTridiagonalSolver::sub_diagonal(const int index) const { - assert(index >= 0); - assert(index < this->matrix_dimension_ - 1); - return this->sub_diagonal_values_[index]; + assert(0 <= index && index < matrix_dimension_ - 1); + return sub_diagonal_values_[index]; } template T& SymmetricTridiagonalSolver::sub_diagonal(const int index) { - assert(index >= 0); - assert(index < this->matrix_dimension_ - 1); - return this->sub_diagonal_values_[index]; + assert(0 <= index && index < matrix_dimension_ - 1); + return sub_diagonal_values_[index]; } template const T& SymmetricTridiagonalSolver::cyclic_corner_element() const { assert(is_cyclic_); - return this->cyclic_corner_element_; + return cyclic_corner_element_; } template T& SymmetricTridiagonalSolver::cyclic_corner_element() { assert(is_cyclic_); - return this->cyclic_corner_element_; + return cyclic_corner_element_; } template @@ -331,18 +255,6 @@ void SymmetricTridiagonalSolver::solveInPlace(T* sol_rhs, T* temp1, T* temp2) template void SymmetricTridiagonalSolver::solveSymmetricTridiagonal(T* x, T* scratch) { - /* ---------------------------------------------------------- */ - /* Based on Cholesky Decomposition: A = L * D * L^T - * - * This function performs Cholesky decomposition on a - * symmetric tridiagonal matrix, factorizing it into - * a lower triangular matrix (L) and a diagonal matrix (D). - * - * By storing the decomposition, this approach enhances - * efficiency for repeated solutions, as matrix factorizations - * need not be recalculated each time. - * ---------------------------------------------------------- */ - // Cholesky Decomposition if (!factorized_) { for (int i = 1; i < matrix_dimension_; i++) { @@ -365,64 +277,15 @@ void SymmetricTridiagonalSolver::solveSymmetricTridiagonal(T* x, T* scratch) for (int i = matrix_dimension_ - 2; i >= 0; i--) { x[i] -= sub_diagonal(i) * x[i + 1]; } - - /* --------------------------------------------------------------- */ - /* Thomas Algorithm: An alternative approach for solving - * tridiagonal systems that does not overwrite the matrix data. - * The algorithm is more stable than the previous approach. - * We enhances numerical stability by performing division directly - * rather than multiplying by the inverse. - * This reduces the risk of precision errors during computation. - * Requires additional 'scratch' storage. - * --------------------------------------------------------------- */ - - // assert(!equals(main_diagonal(0), 0.0)); - // scratch[0] = sub_diagonal(0) / main_diagonal(0); - // x[0] /= main_diagonal(0); - - // for (int i = 1; i < matrix_dimension_ - 1; i++) - // { - // const double divisor = main_diagonal(i) - sub_diagonal(i - 1) * scratch[i - 1]; - // assert(!equals(divisor, 0.0)); - // scratch[i] = sub_diagonal(i) / divisor; - // x[i] = (x[i] - sub_diagonal(i - 1) * x[i - 1]) / divisor; - // } - - // const int i = matrix_dimension_ - 1; - // const double divisor = main_diagonal(i) - sub_diagonal(i - 1) * scratch[i - 1]; - // assert(!equals(divisor, 0.0)); - // x[i] = (x[i] - sub_diagonal(i - 1) * x[i - 1]) / divisor; - - // for (int i = matrix_dimension_ - 2; i >= 0; i--) - // { - // x[i] -= scratch[i] * x[i + 1]; - // } } // ------------------------- // // Cyclic Tridiagonal Solver // // ------------------------- // -/* - * This algorithm implements the Tridiagonal Matrix Algorithm (TDMA) for solving - * symmetric tridiagonal systems of equations, specifically designed to handle - * cyclic boundary conditions. The implementation is based on principles outlined - * in the following resource: - * https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm. - */ - template void SymmetricTridiagonalSolver::solveSymmetricCyclicTridiagonal(T* x, T* u, T* scratch) { - /* ---------------------------------------------------------- */ - /* Cholesky Decomposition: A = L * D * L^T - * This step factorizes the tridiagonal matrix into lower - * triangular (L) and diagonal (D) matrices. While this - * approach may be slightly less stable, it can offer improved - * performance for repeated solves due to the factorization - * being stored internally. - * ---------------------------------------------------------- */ - // Cholesky Decomposition if (!factorized_) { // Shermann-Morrison Adjustment @@ -463,50 +326,4 @@ void SymmetricTridiagonalSolver::solveSymmetricCyclicTridiagonal(T* x, T* u, for (int i = 0; i < matrix_dimension_; i++) { x[i] -= factor * u[i]; } - - /* --------------------------------------------------------------- */ - /* Thomas Algorithm: An alternative approach for solving - * tridiagonal systems that does not overwrite the matrix data. - * The algorithm is more stable than the previous approach. - * We enhances numerical stability by performing division directly - * rather than multiplying by the inverse. - * This reduces the risk of precision errors during computation. - * Requires additional storage in scratch. - * --------------------------------------------------------------- */ - - // const double gamma = -main_diagonal(0); - // const double first_main_diagonal = main_diagonal(0) - gamma; - // const double last_main_diagonal = main_diagonal(matrix_dimension_ - 1) - cyclic_corner_element() * cyclic_corner_element() / gamma; - - // scratch[0] = sub_diagonal(0) / first_main_diagonal; - // x[0] /= first_main_diagonal; - // u[0] = gamma / first_main_diagonal; - - // for (int i = 1; i < matrix_dimension_ - 1; i++) - // { - // const double divisor = main_diagonal(i) - sub_diagonal(i - 1) * scratch[i - 1]; - // scratch[i] = sub_diagonal(i) / divisor; - // x[i] = (x[i] - sub_diagonal(i - 1) * x[i - 1]) / divisor; - // u[i] = (0.0 - sub_diagonal(i - 1) * u[i - 1]) / divisor; - // } - - // const int i = matrix_dimension_ - 1; - // const double divisor = last_main_diagonal - sub_diagonal(i - 1) * scratch[i - 1]; - // x[i] = (x[i] - sub_diagonal(i - 1) * x[i - 1]) / divisor; - // u[i] = (cyclic_corner_element() - sub_diagonal(i - 1) * u[i - 1]) / divisor; - - // for (int i = matrix_dimension_ - 2; i >= 0; i--) - // { - // x[i] -= scratch[i] * x[i + 1]; - // u[i] -= scratch[i] * u[i + 1]; - // } - - // const double dot_product_x_v = x[0] + cyclic_corner_element() / gamma * x[matrix_dimension_ - 1]; - // const double dot_product_u_v = u[0] + cyclic_corner_element() / gamma * u[matrix_dimension_ - 1]; - // const double factor = dot_product_x_v / (1.0 + dot_product_u_v); - - // for (int i = 0; i < matrix_dimension_; i++) - // { - // x[i] -= factor * u[i]; - // } } diff --git a/tests/LinearAlgebra/cyclic_tridiagonal_solver.cpp b/tests/LinearAlgebra/cyclic_tridiagonal_solver.cpp index 568eb0ef..1b8384a6 100644 --- a/tests/LinearAlgebra/cyclic_tridiagonal_solver.cpp +++ b/tests/LinearAlgebra/cyclic_tridiagonal_solver.cpp @@ -93,18 +93,20 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -118,14 +120,14 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_100) @@ -143,18 +145,20 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_100) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -168,14 +172,14 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_100) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_1000) @@ -193,18 +197,20 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_1000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -218,14 +224,14 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_1000) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10000) @@ -243,18 +249,20 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -268,14 +276,14 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10000) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LOW_PRECISION_n_10000) @@ -295,18 +303,20 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagona std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -320,14 +330,14 @@ TEST(ZeroCyclicSymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagona Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } // ------------------------------ // @@ -417,18 +427,21 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - solver.cyclic_corner_element() = dis(gen); - - const SymmetricTridiagonalSolver copy_solver = solver; + solver.cyclic_corner_element() = dis(gen); + const double original_cyclic_corner_element = solver.cyclic_corner_element(); // Fill rhs with random values Vector rhs("rhs", n); @@ -443,16 +456,16 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1] + - copy_solver.cyclic_corner_element() * rhs[n - 1], + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1] + + original_cyclic_corner_element * rhs[n - 1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1] + - copy_solver.cyclic_corner_element() * rhs[0], + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1] + + original_cyclic_corner_element * rhs[0], copy_rhs[n - 1], precision); } @@ -470,18 +483,21 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_100) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - solver.cyclic_corner_element() = dis(gen); - - const SymmetricTridiagonalSolver copy_solver = solver; + solver.cyclic_corner_element() = dis(gen); + const double original_cyclic_corner_element = solver.cyclic_corner_element(); // Fill rhs with random values Vector rhs("rhs", n); @@ -496,16 +512,16 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_100) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1] + - copy_solver.cyclic_corner_element() * rhs[n - 1], + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1] + + original_cyclic_corner_element * rhs[n - 1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1] + - copy_solver.cyclic_corner_element() * rhs[0], + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1] + + original_cyclic_corner_element * rhs[0], copy_rhs[n - 1], precision); } @@ -523,18 +539,21 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_1000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - solver.cyclic_corner_element() = dis(gen); - - const SymmetricTridiagonalSolver copy_solver = solver; + solver.cyclic_corner_element() = dis(gen); + const double original_cyclic_corner_element = solver.cyclic_corner_element(); // Fill rhs with random values Vector rhs("rhs", n); @@ -549,16 +568,16 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_1000) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1] + - copy_solver.cyclic_corner_element() * rhs[n - 1], + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1] + + original_cyclic_corner_element * rhs[n - 1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1] + - copy_solver.cyclic_corner_element() * rhs[0], + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1] + + original_cyclic_corner_element * rhs[0], copy_rhs[n - 1], precision); } @@ -576,18 +595,21 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - solver.cyclic_corner_element() = dis(gen); - - const SymmetricTridiagonalSolver copy_solver = solver; + solver.cyclic_corner_element() = dis(gen); + const double original_cyclic_corner_element = solver.cyclic_corner_element(); // Fill rhs with random values Vector rhs("rhs", n); @@ -602,16 +624,16 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_n_10000) Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1] + - copy_solver.cyclic_corner_element() * rhs[n - 1], + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1] + + original_cyclic_corner_element * rhs[n - 1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1] + - copy_solver.cyclic_corner_element() * rhs[0], + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1] + + original_cyclic_corner_element * rhs[0], copy_rhs[n - 1], precision); } @@ -631,18 +653,21 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LO std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - solver.cyclic_corner_element() = subdiagonal_boost * dis(gen); - - const SymmetricTridiagonalSolver copy_solver = solver; + solver.cyclic_corner_element() = subdiagonal_boost * dis(gen); + const double original_cyclic_corner_element = solver.cyclic_corner_element(); // Fill rhs with random values Vector rhs("rhs", n); @@ -657,15 +682,15 @@ TEST(CyclicSymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LO Vector temp2("temp2", n); solver.solveInPlace(rhs.data(), temp1.data(), temp2.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1] + - copy_solver.cyclic_corner_element() * rhs[n - 1], + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1] + + original_cyclic_corner_element * rhs[n - 1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1] + - copy_solver.cyclic_corner_element() * rhs[0], + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1] + + original_cyclic_corner_element * rhs[0], copy_rhs[n - 1], precision); -} +} \ No newline at end of file diff --git a/tests/LinearAlgebra/tridiagonal_solver.cpp b/tests/LinearAlgebra/tridiagonal_solver.cpp index bee806f4..6d99d87c 100644 --- a/tests/LinearAlgebra/tridiagonal_solver.cpp +++ b/tests/LinearAlgebra/tridiagonal_solver.cpp @@ -130,7 +130,6 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10) { const int n = 10; const double precision = 1e-9; - SymmetricTridiagonalSolver solver(n); solver.is_cyclic(false); @@ -139,18 +138,20 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -163,14 +164,14 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10) Vector temp("temp", n); solver.solveInPlace(rhs.data(), temp.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_100) @@ -187,18 +188,20 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_100) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -211,14 +214,14 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_100) Vector temp("temp", n); solver.solveInPlace(rhs.data(), temp.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_1000) @@ -235,18 +238,20 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_1000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -259,14 +264,14 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_1000) Vector temp("temp", n); solver.solveInPlace(rhs.data(), temp.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10000) @@ -283,18 +288,20 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10000) std::mt19937 gen(42); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = dis(gen); + solver.sub_diagonal(i) = dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -307,24 +314,21 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_n_10000) Vector temp("temp", n); solver.solveInPlace(rhs.data(), temp.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); } TEST(SymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LOW_PRECISION_n_10000) { - const int n = 10000; - - const double precision = 1e-3; - + const int n = 10000; + const double precision = 1e-3; const double subdiagonal_boost = 10000.0; - SymmetricTridiagonalSolver solver(n); solver.is_cyclic(false); @@ -333,18 +337,20 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LOW_PREC std::mt19937 gen(79); std::uniform_real_distribution<> dis(-100.0, 100.0); - // Fill main_diagonal with random values + // Fill main_diagonal with random values and store copies + Vector original_main_diagonal("original_main_diagonal", n); for (int i = 0; i < n; ++i) { - solver.main_diagonal(i) = dis(gen); + solver.main_diagonal(i) = dis(gen); + original_main_diagonal[i] = solver.main_diagonal(i); } - // Fill sub_diagonal with random values + // Fill sub_diagonal with random values and store copies + Vector original_sub_diagonal("original_sub_diagonal", n - 1); for (int i = 0; i < n - 1; ++i) { - solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + solver.sub_diagonal(i) = subdiagonal_boost * dis(gen); + original_sub_diagonal[i] = solver.sub_diagonal(i); } - const SymmetricTridiagonalSolver copy_solver = solver; - // Fill rhs with random values Vector rhs("rhs", n); for (int i = 0; i < n; ++i) { @@ -357,12 +363,12 @@ TEST(SymmetricTridiagonalSolver, random_tridiagonal_boosted_subdiagonal_LOW_PREC Vector temp("temp", n); solver.solveInPlace(rhs.data(), temp.data()); - EXPECT_NEAR(copy_solver.main_diagonal(0) * rhs[0] + copy_solver.sub_diagonal(0) * rhs[1], copy_rhs[0], precision); + EXPECT_NEAR(original_main_diagonal[0] * rhs[0] + original_sub_diagonal[0] * rhs[1], copy_rhs[0], precision); for (int i = 1; i < n - 1; ++i) { - EXPECT_NEAR(copy_solver.sub_diagonal(i - 1) * rhs[i - 1] + copy_solver.main_diagonal(i) * rhs[i] + - copy_solver.sub_diagonal(i) * rhs[i + 1], + EXPECT_NEAR(original_sub_diagonal[i - 1] * rhs[i - 1] + original_main_diagonal[i] * rhs[i] + + original_sub_diagonal[i] * rhs[i + 1], copy_rhs[i], precision); } - EXPECT_NEAR(copy_solver.sub_diagonal(n - 2) * rhs[n - 2] + copy_solver.main_diagonal(n - 1) * rhs[n - 1], - copy_rhs[n - 1], precision); -} + EXPECT_NEAR(original_sub_diagonal[n - 2] * rhs[n - 2] + original_main_diagonal[n - 1] * rhs[n - 1], copy_rhs[n - 1], + precision); +} \ No newline at end of file