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 =
new Matrix(num_rows, num_cols,
true);
110 d_W =
new Matrix(num_rows, num_cols,
true);
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)),
268 CAROM_ASSERT(
d_S != 0);
275 std::cout <<
"Getting snapshot matrix for incremental SVD not implemented" <<
282 double* u,
bool add_without_increase)
284 CAROM_VERIFY(u != 0);
301 if(
d_rank == 0) printf(
"linearly dependent sample!\n");
310 bool linearly_dependent_sample;
313 std::cout <<
"linearly dependent sample! k = " << k <<
"\n";
317 linearly_dependent_sample =
true;
320 linearly_dependent_sample =
true;
329 linearly_dependent_sample =
true;
332 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) {
367 for (
int i = 0; i <
d_dim; ++i) {
398 CAROM_VERIFY(l != 0);
419 Q[q_idx] = l->
item(row);
436 CAROM_VERIFY(A != 0);
444 S->
item(row, col) = 0.0;
457 int lwork = m*(4*m + 7);
458 double* work =
new double [lwork];
481 S->
item(i, i) = sigma[i];
488 double tmp = U->
item(row, col);
489 U->
item(row, col) = U->
item(col, row);
490 U->
item(col, row) = tmp;
514 CAROM_ASSERT(m != 0);
521 for (
int i = 0; i < num_rows; ++i) {
522 tmp += m->
item(i, 0) * m->
item(i, last_col);
525 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.
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...
double checkOrthogonality(const Matrix *m)
Computes and returns the orthogonality of m.
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.
virtual void addNewSample(const Vector *j, const Matrix *A, const Matrix *W, Matrix *sigma)=0
Add a new, unique sample to the SVD.
static const int COMMUNICATE_U
MPI message tag.
virtual void computeBasis()=0
Computes the current basis vectors.
virtual const Vector * getSingularValues()
Returns the singular values for the current time interval.
int d_max_basis_dimension
The maximum basis dimension.
long int d_total_dim
The total dimension of the system.
std::string d_state_file_name
The name of file to which state is saved or restored from.
virtual const Matrix * getTemporalBasis()
Returns the temporal basis vectors for the current time interval as a Matrix.
void constructQ(double *&Q, const Vector *l, double k)
Construct the matrix Q whose SVD is needed.
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.
int d_size
The total number of processors.
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.
virtual const Matrix * getSnapshotMatrix()
Returns the snapshot matrix for the current time interval.
virtual const Matrix * getSpatialBasis()
Returns the basis vectors for the current time interval as a Matrix.
virtual void addLinearlyDependentSample(const Matrix *A, const Matrix *W, const Matrix *sigma)=0
Add a linearly dependent sample to the SVD.
int d_rank
The rank of the processor owning this object.
bool distributed() const
Returns true if the Matrix is distributed.
int numColumns() const
Returns the number of columns in the Matrix. This method will return the same value from each process...
Matrix * transposeMult(const Matrix &other) const
Multiplies the transpose of this Matrix with other and returns the product, reference version.
Matrix * mult(const Matrix &other) const
Multiplies this Matrix with other and returns the product, reference version.
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.
Matrix * d_basis
The globalized basis vectors for the current time interval.
Matrix * d_snapshots
The globalized snapshot vectors for the current time interval.
Matrix * d_basis_right
The globalized right basis vectors for the current time interval.
bool isFirstSample() const
Returns true if the next sample will result in a new time interval.
Vector * d_S
The vector S which is small.
int d_num_samples
Number of samples stored for the current time interval.
Matrix * d_U
The matrix U which is large.
bool d_debug_algorithm
Flag to indicate if results of algorithm should be printed for debugging purposes.
const int d_dim
Dimension of the system.
Matrix * d_W
The matrix U which is large.
double norm() const
Form the norm of this.
Vector * minus(const Vector &other) const
Subtracts other and this and returns the result, reference version.
const double & item(int i) const
Const Vector member access.
double inner_product(const Vector &other) const
Inner product, reference form.
int dim() const
Returns the dimension of the Vector on this processor.