14 #include "IncrementalSVD.h"
15 #include "utils/HDFDatabase.h"
26 #define dgesdd CAROM_FC_GLOBAL(dgesdd, DGESDD)
29 void dgesdd(
char*,
int*,
int*,
double*,
int*,
30 double*,
double*,
int*,
double*,
int*,
31 double*,
int*,
int*,
int*);
38 IncrementalSVD::IncrementalSVD(
40 const std::string& basis_file_name) :
42 d_linearity_tol(options.linearity_tol),
43 d_skip_linearly_dependent(options.skip_linearly_dependent),
44 d_max_basis_dimension(options.max_basis_dimension),
46 d_save_state(options.save_state),
47 d_update_right_SV(options.update_right_SV),
56 MPI_Initialized(&mpi_init);
58 MPI_Comm_size(MPI_COMM_WORLD, &
d_size);
59 MPI_Comm_rank(MPI_COMM_WORLD, &
d_rank);
79 for (
int i = 0; i <
d_size; ++i) {
86 std::ostringstream tmp;
87 tmp << basis_file_name <<
".state." <<
88 std::setw(6) << std::setfill(
'0') <<
d_rank;
101 d_U.reset(
new Matrix(num_rows, num_cols,
true));
110 d_W.reset(
new Matrix(num_rows, num_cols,
true));
143 int num_rows =
d_U->numRows();
145 int num_cols =
d_U->numColumns();
150 int num_dim =
d_S->dim();
158 num_rows =
d_W->numRows();
160 num_cols =
d_W->numColumns();
176 bool add_without_increase)
178 CAROM_VERIFY(u_in != 0);
182 if (u_vec.
norm() == 0.0) {
201 printf(
"%.16e ",
d_S->item(col));
207 for (
int row = 0; row <
d_dim; ++row) {
209 printf(
"%.16e ", basis->item(row, col));
215 for (
int proc = 1; proc <
d_size; ++proc) {
226 for (
int row = 0; row <
d_proc_dims[proc]; ++row) {
228 printf(
"%.16e ", m[idx++]);
234 printf(
"============================================================\n");
239 MPI_Isend(
const_cast<double*
>(&basis->item(0, 0)),
251 std::shared_ptr<const Matrix>
258 std::shared_ptr<const Matrix>
265 std::shared_ptr<const Vector>
268 CAROM_ASSERT(
d_S != 0);
272 std::shared_ptr<const Matrix>
275 std::cout <<
"Getting snapshot matrix for incremental SVD not implemented" <<
282 double* u,
bool add_without_increase)
284 CAROM_VERIFY(u != 0);
289 d_basis->transposeMult(u_vec, l);
298 std::unique_ptr<Vector> e_proj = u_vec.
minus(basisl);
299 double k = e_proj->inner_product(*e_proj);
302 if(
d_rank == 0) printf(
"linearly dependent sample!\n");
311 bool linearly_dependent_sample;
314 std::cout <<
"linearly dependent sample! k = " << k <<
"\n";
318 linearly_dependent_sample =
true;
321 linearly_dependent_sample =
true;
330 linearly_dependent_sample =
true;
333 linearly_dependent_sample =
false;
344 bool result =
svd(Q, A, sigma, W);
358 if(
d_rank == 0) std::cout <<
"adding linearly dependent sample!\n";
362 else if (!linearly_dependent_sample) {
366 std::unique_ptr<Vector> j = u_vec.
minus(basisl);
367 for (
int i = 0; i <
d_dim; ++i) {
407 Q[q_idx] =
d_S->item(col);
415 Q[q_idx] = l.
item(row);
432 CAROM_VERIFY(A != 0);
440 S->
item(row, col) = 0.0;
453 int lwork = m*(4*m + 7);
454 double* work =
new double [lwork];
477 S->
item(i, i) = sigma[i];
484 double tmp = U->
item(row, col);
485 U->
item(row, col) = U->
item(col, row);
486 U->
item(col, row) = tmp;
515 for (
int i = 0; i < num_rows; ++i) {
516 tmp += m.
item(i, 0) * m.
item(i, last_col);
519 MPI_Allreduce(&tmp, &result, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
void getInteger(const std::string &key, int &data)
Reads an integer associated with the supplied key from the currently open database file.
virtual bool close()=0
Closes the currently open database file.
virtual void getDoubleArray(const std::string &key, double *data, int nelements, const bool distributed=false)=0
Reads an array of doubles associated with the supplied key from the currently open database file.
void putInteger(const std::string &key, int data)
Writes an integer associated with the supplied key to currently open database file.
virtual void putDoubleArray(const std::string &key, const double *const data, int nelements, const bool distributed=false)=0
Writes an array of doubles associated with the supplied key to the currently open database file.
virtual bool open(const std::string &file_name, const std::string &type, const MPI_Comm comm=MPI_COMM_NULL)
Opens an existing database file with the supplied name.
std::shared_ptr< const Matrix > getSpatialBasis() override
Returns the basis vectors for the current time interval as a Matrix.
bool d_save_state
If true the state of the SVD will be written to disk when the object is deleted. If there are multipl...
std::shared_ptr< const Matrix > getTemporalBasis() override
Returns the temporal basis vectors for the current time interval as a Matrix.
std::vector< int > d_proc_dims
The dimension of the system on each processor.
Database * d_state_database
Pointer to the database that will hold saved state data if the state is to be saved.
bool d_update_right_SV
Whether to update the right singular vectors.
double d_linearity_tol
Tolerance to determine whether or not a sample is linearly dependent.
static const int COMMUNICATE_U
MPI message tag.
virtual void computeBasis()=0
Computes the current basis vectors.
int d_max_basis_dimension
The maximum basis dimension.
long int d_total_dim
The total dimension of the system.
std::shared_ptr< const Vector > getSingularValues() override
Returns the singular values for the current time interval.
virtual void addLinearlyDependentSample(const Matrix &A, const Matrix &W, const Matrix &sigma)=0
Add a linearly dependent sample to the SVD.
std::string d_state_file_name
The name of file to which state is saved or restored from.
std::shared_ptr< const Matrix > getSnapshotMatrix() override
Returns the snapshot matrix for the current time interval.
virtual bool buildIncrementalSVD(double *u, bool add_without_increase=false)
Adds the new sampled the state vector, u, to the system.
int numSamples()
The number of samples stored.
virtual void addNewSample(const Vector &j, const Matrix &A, const Matrix &W, const Matrix &sigma)=0
Add a new, unique sample to the SVD.
int d_size
The total number of processors.
void constructQ(double *&Q, const Vector &l, double k)
Construct the matrix Q whose SVD is needed.
bool svd(double *A, Matrix *&U, Matrix *&S, Matrix *&V)
Given a matrix, A, returns 2 of the 3 components of its singular value decomposition....
virtual bool takeSample(double *u_in, bool add_without_increase=false)
Sample new state, u_in, at the given time.
virtual ~IncrementalSVD()
Destructor.
virtual void buildInitialSVD(double *u)=0
Constructs the first SVD.
bool d_skip_linearly_dependent
Whether to skip linearly dependent samples.
int d_rank
The rank of the processor owning this object.
double checkOrthogonality(const Matrix &m)
Computes and returns the orthogonality of m.
bool distributed() const
Returns true if the Matrix is distributed.
const double & item(int row, int col) const
Const Matrix member access. Matrix data is stored in row-major format.
int numRows() const
Returns the number of rows of the Matrix on this processor.
bool restore_state
If true the state of the incremental SVD will be restored when the object is created.
double linearity_tol
The tolerance of the incremental SVD algorithm to determine whether or not a sample is linearly depen...
bool save_state
If true the state of the incremental SVD will be written to disk when the object is deleted....
int max_basis_dimension
The maximum dimension of the basis.
std::shared_ptr< Matrix > d_basis_right
The globalized right basis vectors for the current time interval.
std::shared_ptr< Matrix > d_basis
The globalized basis vectors for the current time interval.
bool isFirstSample() const
Returns true if the next sample will result in a new time interval.
int d_num_samples
Number of samples stored for the current time interval.
bool d_debug_algorithm
Flag to indicate if results of algorithm should be printed for debugging purposes.
std::shared_ptr< Matrix > d_snapshots
The globalized snapshot vectors for the current time interval.
std::shared_ptr< Matrix > d_W
The matrix U which is large.
const int d_dim
Dimension of the system.
std::shared_ptr< Vector > d_S
The vector S which is small.
std::shared_ptr< Matrix > d_U
The matrix U which is large.
double norm() const
Form the norm of this.
std::unique_ptr< Vector > minus(const Vector &other) const
Subtracts other and this and returns the result.
const double & item(int i) const
Const Vector member access.
int dim() const
Returns the dimension of the Vector on this processor.