From d3ef2cd2d0db8f9a06b1169aaa2c4a3e11be88a0 Mon Sep 17 00:00:00 2001 From: SeungWon Suh <61690747+swsuh28@users.noreply.github.com> Date: Fri, 21 Jul 2023 00:30:49 -0700 Subject: [PATCH] IncrementalSVDFastUpdate with faster updates + bug fix in IncrementalSVD (#227) * Update IncrementalSVD for faster updates * created a separate class for new features * minor fix in CMakeLists.txt * another minor fix in CMakeLists * add unit test for IncrementalSVDBrand * fixed mpi bug in test_IncrementalSVDBrand * minor styling fix * update args order in unit test --------- Co-authored-by: swsuh28 --- .github/workflows/run_tests/action.yml | 3 + CMakeLists.txt | 1 + lib/CMakeLists.txt | 1 + lib/linalg/BasisGenerator.cpp | 9 +- lib/linalg/Options.h | 8 + lib/linalg/svd/IncrementalSVD.cpp | 13 +- lib/linalg/svd/IncrementalSVDBrand.cpp | 523 ++++++++++++++++++++++++ lib/linalg/svd/IncrementalSVDBrand.h | 188 +++++++++ unit_tests/test_IncrementalSVDBrand.cpp | 147 +++++++ 9 files changed, 886 insertions(+), 7 deletions(-) create mode 100644 lib/linalg/svd/IncrementalSVDBrand.cpp create mode 100644 lib/linalg/svd/IncrementalSVDBrand.h create mode 100644 unit_tests/test_IncrementalSVDBrand.cpp diff --git a/.github/workflows/run_tests/action.yml b/.github/workflows/run_tests/action.yml index d790ed811..34471de9b 100644 --- a/.github/workflows/run_tests/action.yml +++ b/.github/workflows/run_tests/action.yml @@ -22,6 +22,9 @@ runs: mpirun -n 3 --oversubscribe tests/test_RandomizedSVD ./tests/test_StaticSVD mpirun -n 3 --oversubscribe tests/test_StaticSVD + ./tests/test_IncrementalSVDBrand + mpirun -n 3 --oversubscribe tests/test_IncrementalSVDBrand + shell: bash - name: Run regression tests diff --git a/CMakeLists.txt b/CMakeLists.txt index 95d019196..af9036360 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -264,6 +264,7 @@ if(GTEST_FOUND) StaticSVD RandomizedSVD IncrementalSVD + IncrementalSVDBrand GreedyCustomSampler) foreach(stem IN LISTS unit_test_stems) add_executable(test_${stem} unit_tests/test_${stem}.cpp) diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 225081a85..b1128a273 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -30,6 +30,7 @@ set(module_list linalg/svd/IncrementalSVD linalg/svd/IncrementalSVDFastUpdate linalg/svd/IncrementalSVDStandard + linalg/svd/IncrementalSVDBrand linalg/svd/RandomizedSVD linalg/svd/SVD linalg/svd/StaticSVD diff --git a/lib/linalg/BasisGenerator.cpp b/lib/linalg/BasisGenerator.cpp index 1b1020f43..bfd46bbb7 100644 --- a/lib/linalg/BasisGenerator.cpp +++ b/lib/linalg/BasisGenerator.cpp @@ -19,6 +19,7 @@ #include "svd/RandomizedSVD.h" #include "svd/IncrementalSVDStandard.h" #include "svd/IncrementalSVDFastUpdate.h" +#include "svd/IncrementalSVDBrand.h" namespace CAROM { @@ -65,7 +66,13 @@ BasisGenerator::BasisGenerator( d_dt = options.initial_dt; d_next_sample_time = 0.0; - if (options.fast_update) { + if (options.fast_update_brand) { + d_svd.reset( + new IncrementalSVDBrand( + options, + basis_file_name)); + } + else if (options.fast_update) { d_svd.reset( new IncrementalSVDFastUpdate( options, diff --git a/lib/linalg/Options.h b/lib/linalg/Options.h index 2bc833556..082d4a2b2 100644 --- a/lib/linalg/Options.h +++ b/lib/linalg/Options.h @@ -146,6 +146,7 @@ class Options double sampling_tol_, double max_time_between_samples_, bool fast_update_ = false, + bool fast_update_brand_ = false, bool skip_linearly_dependent_ = false ) { @@ -154,6 +155,7 @@ class Options sampling_tol = sampling_tol_; max_time_between_samples = max_time_between_samples_; fast_update = fast_update_; + fast_update_brand = fast_update_brand_; skip_linearly_dependent = skip_linearly_dependent_; return *this; } @@ -299,6 +301,12 @@ class Options */ bool fast_update = false; + /** + * @brief If true use the exact Brand's algorithm for the + * incremental SVD. + */ + bool fast_update_brand = false; + /** * @brief If true skip linearly dependent samples of the * incremental SVD algorithm. diff --git a/lib/linalg/svd/IncrementalSVD.cpp b/lib/linalg/svd/IncrementalSVD.cpp index 36f2f187d..f30551e8d 100644 --- a/lib/linalg/svd/IncrementalSVD.cpp +++ b/lib/linalg/svd/IncrementalSVD.cpp @@ -301,12 +301,13 @@ IncrementalSVD::buildIncrementalSVD( // basisl = basis * l Vector* basisl = d_basis->mult(l); - // Compute k = sqrt(u.u - 2.0*l.l + basisl.basisl) which is ||u - - // basisl||_{2}. This is the error in the projection of u into the - // reduced order space and subsequent lifting back to the full - // order space. - double k = u_vec.inner_product(u_vec) - 2.0*l->inner_product(l) + - basisl->inner_product(basisl); + // Computing as k = sqrt(u.u - 2.0*l.l + basisl.basisl) + // results in catastrophic cancellation, and must be avoided. + // Instead we compute as k = sqrt((u-basisl).(u-basisl)). + Vector* e_proj = u_vec.minus(basisl); + double k = e_proj->inner_product(e_proj); + delete e_proj; + if (k <= 0) { if(d_rank == 0) printf("linearly dependent sample!\n"); k = 0; diff --git a/lib/linalg/svd/IncrementalSVDBrand.cpp b/lib/linalg/svd/IncrementalSVDBrand.cpp new file mode 100644 index 000000000..beea05218 --- /dev/null +++ b/lib/linalg/svd/IncrementalSVDBrand.cpp @@ -0,0 +1,523 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: The concrete implementation of the incremental SVD algorithm +// using Matthew Brand's "fast update" method. + +#include "IncrementalSVDBrand.h" +#include "utils/HDFDatabase.h" + +#include "mpi.h" + +#include +#include + +namespace CAROM { + +IncrementalSVDBrand::IncrementalSVDBrand( + Options options, + const std::string& basis_file_name) : + IncrementalSVD( + options, + basis_file_name), + d_Up(0), + d_singular_value_tol(options.singular_value_tol) +{ + CAROM_VERIFY(options.singular_value_tol >= 0); + + // If the state of the SVD is to be restored, do it now. The base class, + // IncrementalSVD, has already opened the database and restored the state + // common to all incremental algorithms. This particular class must also + // read the state of d_Up and then compute the basis. If the database could + // not be found then we can not restore the state. + if (options.restore_state && d_state_database) { + // Read d_Up. + int num_rows; + d_state_database->getInteger("Up_num_rows", num_rows); + int num_cols; + d_state_database->getInteger("Up_num_cols", num_cols); + d_Up = new Matrix(num_rows, num_cols, true); + d_state_database->getDoubleArray("Up", + &d_Up->item(0, 0), + num_rows*num_cols); + + // Close and delete the database. + d_state_database->close(); + delete d_state_database; + + // Compute the basis. + computeBasis(); + } +} + +IncrementalSVDBrand::~IncrementalSVDBrand() +{ + // If the state of the SVD is to be saved, then create the database now. + // The IncrementalSVD base class destructor will save d_S and d_U. This + // derived class must save its specific state data, d_Up. + // + // If there are multiple time intervals then saving and restoring the state + // does not make sense as there is not one, all encompassing, basis. + if (d_save_state && d_time_interval_start_times.size() == 1) { + // Create state database file. + d_state_database = new HDFDatabase(); + d_state_database->create(d_state_file_name); + + // Save d_Up. + int num_rows = d_Up->numRows(); + d_state_database->putInteger("Up_num_rows", num_rows); + int num_cols = d_Up->numColumns(); + d_state_database->putInteger("Up_num_cols", num_cols); + d_state_database->putDoubleArray("Up", + &d_Up->item(0, 0), + num_rows*num_cols); + } + + // Delete data members. + if (d_Up) { + delete d_Up; + } +} + +const Matrix* +IncrementalSVDBrand::getSpatialBasis() +{ + updateSpatialBasis(); // WARNING: this is costly + + CAROM_ASSERT(d_basis != 0); + return d_basis; +} + +const Matrix* +IncrementalSVDBrand::getTemporalBasis() +{ + updateTemporalBasis(); + + CAROM_ASSERT(d_basis_right != 0); + return d_basis_right; +} + +void +IncrementalSVDBrand::buildInitialSVD( + double* u, + double time) +{ + CAROM_VERIFY(u != 0); + CAROM_VERIFY(time >= 0.0); + + // We have a new time interval. + + // If this is not the first time interval then delete d_basis, d_U, d_Up, + // and d_S of the just completed time interval. + int num_time_intervals = + static_cast(d_time_interval_start_times.size()); + if (num_time_intervals > 0) { + delete d_basis; + delete d_U; + delete d_Up; + delete d_S; + delete d_W; + } + increaseTimeInterval(); + d_time_interval_start_times[num_time_intervals] = time; + + // Build d_S for this new time interval. + d_S = new Vector(1, false); + Vector u_vec(u, d_dim, true); + double norm_u = u_vec.norm(); + d_S->item(0) = norm_u; + + // Build d_Up for this new time interval. + d_Up = new Matrix(1, 1, false); + d_Up->item(0, 0) = 1.0; + + // Build d_U for this new time interval. + d_U = new Matrix(d_dim, 1, true); + for (int i = 0; i < d_dim; ++i) { + d_U->item(i, 0) = u[i]/norm_u; + } + + // Build d_W for this new time interval. + if (d_update_right_SV) { + d_W = new Matrix(1, 1, false); + d_W->item(0, 0) = 1.0; + } + + // We now have the first sample for the new time interval. + d_num_samples = 1; + d_num_rows_of_W = 1; + +} + +bool +IncrementalSVDBrand::buildIncrementalSVD( + double* u, bool add_without_increase) +{ + CAROM_VERIFY(u != 0); + + // Compute the projection error + // (accurate down to the machine precision) + Vector u_vec(u, d_dim, true); + Vector e_proj(u, d_dim, true); + e_proj -= *(d_U->mult(d_U->transposeMult(e_proj))); // Gram-Schmidt + e_proj -= *(d_U->mult(d_U->transposeMult(e_proj))); // Re-orthogonalization + + double k = e_proj.inner_product(e_proj); + if (k <= 0) { + if(d_rank == 0) printf("linearly dependent sample!\n"); + k = 0; + } + else { + k = sqrt(k); + } + + // Use k to see if the vector addressed by u is linearly dependent + // on the left singular vectors. + bool linearly_dependent_sample; + if ( k < d_linearity_tol ) { + if(d_rank == 0) { + std::cout << "linearly dependent sample! k = " << k << "\n"; + std::cout << "d_linearity_tol = " << d_linearity_tol << "\n"; + } + k = 0; + linearly_dependent_sample = true; + } else if ( d_num_samples >= d_max_basis_dimension || add_without_increase ) { + k = 0; + linearly_dependent_sample = true; + } + // Check to see if the "number of samples" (in IncrementalSVD and + // its subclasses, d_num_samples appears to be equal to the number + // of columns of the left singular vectors) is greater than or equal + // to the dimension of snapshot vectors. If so, then the vector + // addressed by the pointer u must be linearly dependent on the left + // singular vectors. + else if (d_num_samples >= d_total_dim) { + linearly_dependent_sample = true; + } + else { + linearly_dependent_sample = false; + } + + // Create Q. + double* Q; + Vector* U_mult_u = new Vector(d_U->transposeMult(u_vec)->getData(), + d_num_samples, + false); + Vector* l = d_Up->transposeMult(U_mult_u); + constructQ(Q, l, k); + delete U_mult_u, l; + + // Now get the singular value decomposition of Q. + Matrix* A; + Matrix* W; + Matrix* sigma; + bool result = svd(Q, A, sigma, W); + + // Done with Q. + delete [] Q; + + // If the svd was successful then add the sample. Otherwise clean up and + // return. + if (result) { + + // We need to add the sample if it is not linearly dependent or if it is + // linearly dependent and we are not skipping linearly dependent samples. + if ((linearly_dependent_sample && !d_skip_linearly_dependent) ) { + // This sample is linearly dependent and we are not skipping linearly + // dependent samples. + if(d_rank == 0) std::cout << "adding linearly dependent sample!\n"; + addLinearlyDependentSample(A, W, sigma); + delete sigma; + } + else if (!linearly_dependent_sample) { + // This sample is not linearly dependent. + + // Compute j + Vector* j = new Vector(e_proj.getData(), d_dim, false); + for (int i = 0; i < d_dim; ++i) { + j->item(i) /= k; + } + + // addNewSample will assign sigma to d_S hence it should not be + // deleted upon return. + addNewSample(j, A, W, sigma); + delete j; + } + delete A; + delete W; + } + else { + delete A; + delete W; + delete sigma; + } + return result; +} + +void +IncrementalSVDBrand::updateSpatialBasis() +{ + d_basis = d_U->mult(d_Up); + + // remove the smallest singular value if it is smaller than d_singular_value_tol + if ( (d_singular_value_tol != 0.0) && + (d_S->item(d_num_samples-1) < d_singular_value_tol) && + (d_num_samples != 1) ) { + + if (d_rank == 0) { + std::cout << + "removing a spatial basis corresponding to the small singular value!\n"; + } + + Matrix* d_basis_new = new Matrix(d_dim, d_num_samples-1, + d_basis->distributed()); + for (int row = 0; row < d_dim; ++row) { + for (int col = 0; col < d_num_samples-1; ++col) { + d_basis_new->item(row, col) = d_basis->item(row,col); + } + } + delete d_basis; + d_basis = d_basis_new; + } + + // Reorthogonalize if necessary. + // (not likely to be called anymore but left for safety) + if (fabs(checkOrthogonality(d_basis)) > + std::numeric_limits::epsilon()*static_cast(d_num_samples)) { + d_basis->orthogonalize(); + } + +} + +void +IncrementalSVDBrand::updateTemporalBasis() +{ + delete d_basis_right; + d_basis_right = new Matrix(*d_W); + + // remove the smallest singular value if it is smaller than d_singular_value_tol + if ( (d_singular_value_tol != 0.0) && + (d_S->item(d_num_samples-1) < d_singular_value_tol) && + (d_num_samples != 1) ) { + + if (d_rank == 0) { + std::cout << + "removing a temporal basis corresponding to the small singular value!\n"; + } + + Matrix* d_basis_right_new = new Matrix(d_num_rows_of_W, d_num_samples-1, + d_basis_right->distributed()); + for (int row = 0; row < d_num_rows_of_W; ++row) { + for (int col = 0; col < d_num_samples-1; ++col) { + d_basis_right_new->item(row, col) = d_basis_right->item(row,col); + } + } + delete d_basis_right; + d_basis_right = d_basis_right_new; + } + + // Reorthogonalize if necessary. + // (not likely to be called anymore but left for safety) + if (fabs(checkOrthogonality(d_basis_right)) > + std::numeric_limits::epsilon()*d_num_samples) { + d_basis_right->orthogonalize(); + } + +} + +void +IncrementalSVDBrand::computeBasis() +{ + if(d_rank == 0) { + std::cout << "d_num_samples = " << d_num_samples << "\n"; + std::cout << "d_num_rows_of_W = " << d_num_rows_of_W << "\n"; + std::cout << "d_singular_value_tol = " << d_singular_value_tol << "\n"; + std::cout << "smallest SV = " << d_S->item(d_num_samples-1) << "\n"; + if (d_num_samples > 1) { + std::cout << "next smallest SV = " << d_S->item(d_num_samples-2) << "\n"; + } + } + + updateSpatialBasis(); + if (d_update_right_SV) + { + updateTemporalBasis(); + } + + // remove the smallest singular value if it is smaller than d_singular_value_tol + if ( (d_singular_value_tol != 0.0) && + (d_S->item(d_num_samples-1) < d_singular_value_tol) && + (d_num_samples != 1) ) { + + --d_num_samples; + } +} + +void +IncrementalSVDBrand::addLinearlyDependentSample( + const Matrix* A, + const Matrix* W, + const Matrix* sigma) +{ + CAROM_VERIFY(A != 0); + CAROM_VERIFY(sigma != 0); + + // Chop a row and a column off of A to form Amod. Also form + // d_S by chopping a row and a column off of sigma. + Matrix Amod(d_num_samples, d_num_samples, false); + for (int row = 0; row < d_num_samples; ++row) { + for (int col = 0; col < d_num_samples; ++col) { + Amod.item(row, col) = A->item(row, col); + if (row == col) + { + d_S->item(col) = sigma->item(row, col); + } + } + } + + // Multiply d_Up and Amod and put result into d_Up. + Matrix* Up_times_Amod = d_Up->mult(Amod); + delete d_Up; + d_Up = Up_times_Amod; + + Matrix* new_d_W; + if (d_update_right_SV) { + // The new d_W is the product of the current d_W extended by another row + // and column and W. The only new value in the extended version of d_W + // that is non-zero is the new lower right value and it is 1. We will + // construct this product without explicitly forming the extended version of + // d_W. + new_d_W = new Matrix(d_num_rows_of_W+1, d_num_samples, false); + for (int row = 0; row < d_num_rows_of_W; ++row) { + for (int col = 0; col < d_num_samples; ++col) { + double new_d_W_entry = 0.0; + for (int entry = 0; entry < d_num_samples; ++entry) { + new_d_W_entry += d_W->item(row, entry)*W->item(entry, col); + } + new_d_W->item(row, col) = new_d_W_entry; + } + } + for (int col = 0; col < d_num_samples; ++col) { + new_d_W->item(d_num_rows_of_W, col) = W->item(d_num_samples, col); + } + delete d_W; + d_W = new_d_W; + ++d_num_rows_of_W; + } + +} + +void +IncrementalSVDBrand::addNewSample( + const Vector* j, + const Matrix* A, + const Matrix* W, + Matrix* sigma) +{ + CAROM_VERIFY(j != 0); + CAROM_VERIFY(A != 0); + CAROM_VERIFY(sigma != 0); + + // Add j as a new column of d_U. + Matrix* newU = new Matrix(d_dim, d_num_samples+1, true); + for (int row = 0; row < d_dim; ++row) { + for (int col = 0; col < d_num_samples; ++col) { + newU->item(row, col) = d_U->item(row, col); + } + newU->item(row, d_num_samples) = j->item(row); + } + delete d_U; + d_U = newU; + + Matrix* new_d_W; + if (d_update_right_SV) { + // The new d_W is the product of the current d_W extended by another row + // and column and W. The only new value in the extended version of d_W + // that is non-zero is the new lower right value and it is 1. We will + // construct this product without explicitly forming the extended version of + // d_W. + new_d_W = new Matrix(d_num_rows_of_W+1, d_num_samples+1, false); + for (int row = 0; row < d_num_rows_of_W; ++row) { + for (int col = 0; col < d_num_samples+1; ++col) { + double new_d_W_entry = 0.0; + for (int entry = 0; entry < d_num_samples; ++entry) { + new_d_W_entry += d_W->item(row, entry)*W->item(entry, col); + } + new_d_W->item(row, col) = new_d_W_entry; + } + } + for (int col = 0; col < d_num_samples+1; ++col) { + new_d_W->item(d_num_rows_of_W, col) = W->item(d_num_samples, col); + } + delete d_W; + d_W = new_d_W; + } + + // The new d_Up is the product of the current d_Up extended by another row + // and column and A. The only new value in the extended version of d_Up + // that is non-zero is the new lower right value and it is 1. We will + // construct this product without explicitly forming the extended version of + // d_Up. + Matrix* new_d_Up = new Matrix(d_num_samples+1, d_num_samples+1, false); + for (int row = 0; row < d_num_samples; ++row) { + for (int col = 0; col < d_num_samples+1; ++col) { + double new_d_Up_entry = 0.0; + for (int entry = 0; entry < d_num_samples; ++entry) { + new_d_Up_entry += d_Up->item(row, entry)*A->item(entry, col); + } + new_d_Up->item(row, col) = new_d_Up_entry; + } + } + for (int col = 0; col < d_num_samples+1; ++col) { + new_d_Up->item(d_num_samples, col) = A->item(d_num_samples, col); + } + delete d_Up; + d_Up = new_d_Up; + + // d_S = sigma. + delete d_S; + int num_dim = std::min(sigma->numRows(), sigma->numColumns()); + d_S = new Vector(num_dim, false); + for (int i = 0; i < num_dim; i++) { + d_S->item(i) = sigma->item(i,i); + } + + // We now have another sample. + ++d_num_samples; + ++d_num_rows_of_W; + + // Reorthogonalize if necessary. + long int max_U_dim; + if (d_num_samples > d_total_dim) { + max_U_dim = d_num_samples; + } + else { + max_U_dim = d_total_dim; + } + if (fabs(checkOrthogonality(d_Up)) > + std::numeric_limits::epsilon()*static_cast(max_U_dim)) { + d_Up->orthogonalize(); + } + if (fabs(checkOrthogonality(d_U)) > + std::numeric_limits::epsilon()*static_cast(max_U_dim)) { + d_U->orthogonalize(); // Will not be called, but just in case + } + + if(d_update_right_SV) + { + if (fabs(checkOrthogonality(d_W)) > + std::numeric_limits::epsilon()*d_num_samples) { + d_W->orthogonalize(); + } + } + +} + +} diff --git a/lib/linalg/svd/IncrementalSVDBrand.h b/lib/linalg/svd/IncrementalSVDBrand.h new file mode 100644 index 000000000..af44beaf4 --- /dev/null +++ b/lib/linalg/svd/IncrementalSVDBrand.h @@ -0,0 +1,188 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: The concrete implementation of the incremental SVD algorithm +// using Matthew Brand's (exact) "fast update" method. + +#ifndef included_IncrementalSVDBrand_h +#define included_IncrementalSVDBrand_h + +#include "IncrementalSVD.h" +#include "linalg/Options.h" + +namespace CAROM { + +/** + * Class IncrementalSVDBrand implements Brand's fast update incremental SVD + * algorithm by implementing the pure virtual methods of the IncrementalSVD + * base class. + */ +class IncrementalSVDBrand : public IncrementalSVD +{ +public: + /** + * @brief Destructor. + */ + ~IncrementalSVDBrand(); + + /** + * @brief Returns the basis vectors for the current time interval as a + * Matrix. + * + * @return The basis vectors for the current time interval. + */ + const Matrix* + getSpatialBasis() override; + + /** + * @brief Returns the temporal basis vectors for the current time interval + * as a Matrix. + * + * @return The temporal basis vectors for the current time interval. + */ + const Matrix* + getTemporalBasis() override; + +private: + friend class BasisGenerator; + + /** + * @brief Constructor. + * + * @param[in] options The struct containing the options for this SVD + * implementation. + * @param[in] basis_file_name The base part of the name of the file + * containing the basis vectors. Each process + * will append its process ID to this base + * name. + * @see Options + */ + IncrementalSVDBrand( + Options options, + const std::string& basis_file_name); + + /** + * @brief Unimplemented default constructor. + */ + IncrementalSVDBrand(); + + /** + * @brief Unimplemented copy constructor. + */ + IncrementalSVDBrand( + const IncrementalSVDBrand& other); + + /** + * @brief Unimplemented assignment operator. + */ + IncrementalSVDBrand& + operator = ( + const IncrementalSVDBrand& rhs); + + /** + * @brief Constructs the first SVD. + * + * @pre u != 0 + * @pre time >= 0.0 + * + * @param[in] u The first state. + * @param[in] time The simulation time for the first state. + */ + virtual + void + buildInitialSVD( + double* u, + double time); + + /** + * @brief Adds the new sampled the state vector, u, to the system. + * + * @pre u != 0 + * + * @param[in] u The new state. + * @param[in] add_without_increase If true, addLinearlyDependent is invoked. + * + * @return True if building the incremental SVD was successful. + */ + bool + buildIncrementalSVD( + double* u, bool add_without_increase = false) override; + + /** + * @brief Update the current spatial basis vectors. + */ + void + updateSpatialBasis(); + + /** + * @brief Update the current temporal basis vectors. + */ + void + updateTemporalBasis(); + + /** + * @brief Computes the current basis vectors. + */ + virtual + void + computeBasis(); + + /** + * @brief Add a linearly dependent sample to the SVD. + * + * @pre A != 0 + * @pre sigma != 0 + * + * @param[in] A The left singular vectors. + * @param[in] W The right singular vectors. + * @param[in] sigma The singular values. + */ + void + addLinearlyDependentSample( + const Matrix* A, + const Matrix* W, + const Matrix* sigma); + + /** + * @brief Add a new, unique sample to the SVD. + * + * @pre j != 0 + * @pre A != 0 + * @pre W != 0 + * @pre sigma != 0 + * + * @param[in] j The new column of d_U. + * @param[in] A The left singular vectors. + * @param[in] W The right singular vectors. + * @param[in] sigma The singular values. + */ + void + addNewSample( + const Vector* j, + const Matrix* A, + const Matrix* W, + Matrix* sigma); + + /** + * @brief The matrix U'. U' is not distributed and the entire matrix + * exists on each processor. + */ + Matrix* d_Up; + + /** + * @brief The tolerance value used to remove small singular values. + */ + double d_singular_value_tol; + +}; + +} + +#endif diff --git a/unit_tests/test_IncrementalSVDBrand.cpp b/unit_tests/test_IncrementalSVDBrand.cpp new file mode 100644 index 000000000..eeaddd1d5 --- /dev/null +++ b/unit_tests/test_IncrementalSVDBrand.cpp @@ -0,0 +1,147 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC + * and other libROM project developers. See the top-level COPYRIGHT + * file for details. + * + * SPDX-License-Identifier: (Apache-2.0 OR MIT) + * + *****************************************************************************/ + +// Description: This source file is a test runner that uses the Google Test +// Framework to run unit tests on the CAROM::IncrementalSVD class. + +#include + +#ifdef CAROM_HAS_GTEST +#include +#include +#include "linalg/BasisGenerator.h" + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +TEST(IncrementalSVDBrandTest, Test_IncrementalSVDBrand) +{ + // Get the rank of this process, and the number of processors. + int mpi_init, d_rank, d_num_procs; + MPI_Initialized(&mpi_init); + if (mpi_init == 0) { + MPI_Init(nullptr, nullptr); + } + + MPI_Comm_rank(MPI_COMM_WORLD, &d_rank); + MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs); + + int num_total_rows = 5; + int d_num_rows = num_total_rows / d_num_procs; + if (num_total_rows % d_num_procs > d_rank) { + d_num_rows++; + } + int *row_offset = new int[d_num_procs + 1]; + row_offset[d_num_procs] = num_total_rows; + row_offset[d_rank] = d_num_rows; + + MPI_Allgather(MPI_IN_PLACE, + 1, + MPI_INT, + row_offset, + 1, + MPI_INT, + MPI_COMM_WORLD); + + for (int i = d_num_procs - 1; i >= 0; i--) { + row_offset[i] = row_offset[i + 1] - row_offset[i]; + } + + double* sample1 = new double[5] {0.5377, 1.8339, -2.2588, 0.8622, 0.3188}; + double* sample2 = new double[5] {-1.3077, -0.4336, 0.3426, 3.5784, 2.7694}; + double* sample3 = new double[5] {-1.3499, 3.0349, 0.7254, -0.0631, 0.7147}; + + double* basis_true_ans = new double[15] { + 3.08158946098238906153E-01, -9.49897947980619661301E-02, -4.50691774108525788911E-01, + -1.43697905723455976457E-01, 9.53289043424090820622E-01, 8.77767692937209131898E-02, + -2.23655845793717528158E-02, -2.10628953513210204207E-01, 8.42235962392685943989E-01, + -7.29903965154318323805E-01, -1.90917141788945754488E-01, -2.77280930877637610266E-01, + -5.92561353877168350834E-01, -3.74570084880578441089E-02, 5.40928141934190823137E-02 + }; + + double* basis_right_true_ans = new double[9] { + -1.78651649346571794741E-01, 5.44387957786310106023E-01, -8.19588518467042281834E-01, + -9.49719639253861602768E-01, -3.13100149275943651084E-01, -9.50441422536040881122E-04, + -2.57130696341890396805E-01, 7.78209514167382598870E-01, 5.72951792961765460355E-01 + }; + + double* sv_true_ans = new double[3] { + 4.84486375065219387892E+00, 3.66719976398777269821E+00, 2.69114625366671811335E+00 + }; + + bool fast_update = true; + bool fast_update_brand = true; + CAROM::Options incremental_svd_options = CAROM::Options(d_num_rows, 3, -1, true) + .setMaxBasisDimension(num_total_rows) + .setIncrementalSVD(1e-1, + 1e-1, + 1e-1, + 1e-1, + fast_update, + fast_update_brand, + false); + + CAROM::BasisGenerator sampler( + incremental_svd_options, + true, + "irrelevant.txt"); + sampler.takeSample(&sample1[row_offset[d_rank]], 0, 1e-1); + sampler.takeSample(&sample2[row_offset[d_rank]], 0, 1e-1); + sampler.takeSample(&sample3[row_offset[d_rank]], 0, 1e-1); + + const CAROM::Matrix* d_basis = sampler.getSpatialBasis(); + const CAROM::Matrix* d_basis_right = sampler.getTemporalBasis(); + const CAROM::Vector* sv = sampler.getSingularValues(); + + EXPECT_EQ(d_basis->numRows(), d_num_rows); + EXPECT_EQ(d_basis->numColumns(), 3); + EXPECT_EQ(d_basis_right->numRows(), 3); + EXPECT_EQ(d_basis_right->numColumns(), 3); + EXPECT_EQ(sv->dim(), 3); + + double* d_basis_vals = d_basis->getData(); + double* d_basis_right_vals = d_basis_right->getData(); + + + for (int i = 0; i < d_num_rows * 3; i++) { + EXPECT_NEAR(abs(d_basis_vals[i]), + abs(basis_true_ans[row_offset[d_rank] * 3 + i]), 1e-7); + } + + for (int i = 0; i < 9; i++) { + EXPECT_NEAR(abs(d_basis_right_vals[i]), abs(basis_right_true_ans[i]), 1e-7); + } + + for (int i = 0; i < 3; i++) { + EXPECT_NEAR(sv->item(i), sv_true_ans[i], 1e-7); + } + +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} +#else // #ifndef CAROM_HAS_GTEST +int main() +{ + std::cout << "libROM was compiled without Google Test support, so unit " + << "tests have been disabled. To enable unit tests, compile " + << "libROM with Google Test support." << std::endl; +} +#endif // #endif CAROM_HAS_GTEST