libROM  v1.0
Data-driven physical simulation library
IncrementalSVD.cpp
1 
11 // Description: The abstract incremental SVD algorithm defines algorithm
12 // interface.
13 
14 #include "IncrementalSVD.h"
15 #include "utils/HDFDatabase.h"
16 
17 #include "mpi.h"
18 
19 #include <cmath>
20 #include <iomanip>
21 #include <limits>
22 #include <stdio.h>
23 #include <sstream>
24 
25 /* Use Autotools-detected Fortran name-mangling scheme */
26 #define dgesdd CAROM_FC_GLOBAL(dgesdd, DGESDD)
27 
28 extern "C" {
29  void dgesdd(char*, int*, int*, double*, int*,
30  double*, double*, int*, double*, int*,
31  double*, int*, int*, int*);
32 }
33 
34 namespace CAROM {
35 
36 const int IncrementalSVD::COMMUNICATE_U = 666;
37 
38 IncrementalSVD::IncrementalSVD(
39  Options options,
40  const std::string& basis_file_name) :
41  SVD(options),
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),
45  d_total_dim(0),
46  d_save_state(options.save_state),
47  d_update_right_SV(options.update_right_SV),
48  d_state_database(0)
49 {
50  CAROM_VERIFY(options.linearity_tol > 0.0);
51  CAROM_VERIFY(options.max_basis_dimension > 0);
52 
53  // Get the number of processors, the dimensions for each process, and the
54  // total dimension.
55  int mpi_init;
56  MPI_Initialized(&mpi_init);
57  if (mpi_init) {
58  MPI_Comm_size(MPI_COMM_WORLD, &d_size);
59  MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
60  }
61  else {
62  d_size = 1;
63  d_rank = 0;
64  }
65  d_proc_dims.reserve(d_size);
66  d_proc_dims.resize(d_size);
67  if (mpi_init) {
68  MPI_Allgather(&d_dim,
69  1,
70  MPI_INT,
71  &d_proc_dims[0],
72  1,
73  MPI_INT,
74  MPI_COMM_WORLD);
75  }
76  else {
77  d_proc_dims[0] = d_dim;
78  }
79  for (int i = 0; i < d_size; ++i) {
81  }
82 
83  // If the state of the SVD is to be restored then open the database and
84  // restore the necessary data from the database now.
85  if (options.save_state || options.restore_state) {
86  std::ostringstream tmp;
87  tmp << basis_file_name << ".state." <<
88  std::setw(6) << std::setfill('0') << d_rank;
89  d_state_file_name = tmp.str();
90  }
91  if (options.restore_state) {
92  // Open state database file.
94  bool is_good = d_state_database->open(d_state_file_name, "r");
95  if (is_good) {
96  // Read d_U.
97  int num_rows;
98  d_state_database->getInteger("U_num_rows", num_rows);
99  int num_cols;
100  d_state_database->getInteger("U_num_cols", num_cols);
101  d_U = new Matrix(num_rows, num_cols, true);
103  &d_U->item(0, 0),
104  num_rows*num_cols);
105 
106  if (d_update_right_SV) {
107  // Read d_W.
108  d_state_database->getInteger("W_num_rows", num_rows);
109  d_state_database->getInteger("W_num_cols", num_cols);
110  d_W = new Matrix(num_rows, num_cols, true);
112  &d_W->item(0, 0),
113  num_rows*num_cols);
114  }
115 
116  // Read d_S.
117  int num_dim;
118  d_state_database->getInteger("S_num_dim", num_dim);
119  d_S = new Vector(num_dim, false);
121  &d_S->item(0),
122  num_dim);
123 
124  // Set d_num_samples.
125  d_num_samples = num_cols;
126  }
127  else {
128  delete d_state_database;
129  d_state_database = 0;
130  }
131  }
132 }
133 
135 {
136  // If the state of the SVD is to be saved, then save d_S and d_U now. The
137  // derived class has already created the database.
138  //
139  // If there are multiple time intervals then saving and restoring the state
140  // does not make sense as there is not one, all encompassing, basis.
141  if (d_save_state && (!isFirstSample())) {
142  // Save d_U.
143  int num_rows = d_U->numRows();
144  d_state_database->putInteger("U_num_rows", num_rows);
145  int num_cols = d_U->numColumns();
146  d_state_database->putInteger("U_num_cols", num_cols);
147  d_state_database->putDoubleArray("U", &d_U->item(0, 0), num_rows*num_cols);
148 
149  // Save d_S.
150  int num_dim = d_S->dim();
151  d_state_database->putInteger("S_num_dim", num_dim);
153  &d_S->item(0),
154  num_dim);
155 
156  if (d_update_right_SV) {
157  // Save d_W.
158  num_rows = d_W->numRows();
159  d_state_database->putInteger("W_num_rows", num_rows);
160  num_cols = d_W->numColumns();
161  d_state_database->putInteger("W_num_cols", num_cols);
163  &d_W->item(0, 0),
164  num_rows*num_cols);
165  }
166 
167  // Close state database file and delete database object.
169  delete d_state_database;
170  }
171 }
172 
173 bool
175  double* u_in,
176  bool add_without_increase)
177 {
178  CAROM_VERIFY(u_in != 0);
179 
180  // Check that u_in is not non-zero.
181  Vector u_vec(u_in, d_dim, true);
182  if (u_vec.norm() == 0.0) {
183  return false;
184  }
185 
186  // If this is the first SVD then build it. Otherwise add this sample to the
187  // system.
188  bool result = true;
189  if (isFirstSample()) {
190  buildInitialSVD(u_in);
191  }
192  else {
193  result = buildIncrementalSVD(u_in,add_without_increase);
194  }
195 
196  if (d_debug_algorithm) {
197  const Matrix* basis = getSpatialBasis();
198  if (d_rank == 0) {
199  // Print d_S.
200  for (int col = 0; col < d_num_samples; ++col) {
201  printf("%.16e ", d_S->item(col));
202  printf("\n");
203  }
204  printf("\n");
205 
206  // Print process 0's part of the basis.
207  for (int row = 0; row < d_dim; ++row) {
208  for (int col = 0; col < d_num_samples; ++col) {
209  printf("%.16e ", basis->item(row, col));
210  }
211  printf("\n");
212  }
213 
214  // Gather other processor's parts of the basis and print them.
215  for (int proc = 1; proc < d_size; ++proc) {
216  double* m = new double[d_proc_dims[proc]*d_num_samples];
217  MPI_Status status;
218  MPI_Recv(m,
220  MPI_DOUBLE,
221  proc,
223  MPI_COMM_WORLD,
224  &status);
225  int idx = 0;
226  for (int row = 0; row < d_proc_dims[proc]; ++row) {
227  for (int col = 0; col < d_num_samples; ++col) {
228  printf("%.16e ", m[idx++]);
229  }
230  printf("\n");
231  }
232  delete [] m;
233  }
234  printf("============================================================\n");
235  }
236  else {
237  // Send this processor's part of the basis to process 0.
238  MPI_Request request;
239  MPI_Isend(const_cast<double*>(&basis->item(0, 0)),
241  MPI_DOUBLE,
242  0,
244  MPI_COMM_WORLD,
245  &request);
246  }
247  }
248  return result;
249 }
250 
251 const Matrix*
253 {
254  CAROM_ASSERT(d_basis != 0);
255  return d_basis;
256 }
257 
258 const Matrix*
260 {
261  CAROM_ASSERT(d_basis != 0);
262  return d_basis_right;
263 }
264 
265 const Vector*
267 {
268  CAROM_ASSERT(d_S != 0);
269  return d_S;
270 }
271 
272 const Matrix*
274 {
275  std::cout << "Getting snapshot matrix for incremental SVD not implemented" <<
276  std::endl;
277  return d_snapshots;
278 }
279 
280 bool
282  double* u, bool add_without_increase)
283 {
284  CAROM_VERIFY(u != 0);
285 
286  // l = basis' * u
287  Vector u_vec(u, d_dim, true);
288  Vector* l = d_basis->transposeMult(u_vec);
289 
290  // basisl = basis * l
291  Vector* basisl = d_basis->mult(l);
292 
293  // Computing as k = sqrt(u.u - 2.0*l.l + basisl.basisl)
294  // results in catastrophic cancellation, and must be avoided.
295  // Instead we compute as k = sqrt((u-basisl).(u-basisl)).
296  Vector* e_proj = u_vec.minus(basisl);
297  double k = e_proj->inner_product(e_proj);
298  delete e_proj;
299 
300  if (k <= 0) {
301  if(d_rank == 0) printf("linearly dependent sample!\n");
302  k = 0;
303  }
304  else {
305  k = sqrt(k);
306  }
307 
308  // Use k to see if the vector addressed by u is linearly dependent
309  // on the left singular vectors.
310  bool linearly_dependent_sample;
311  if ( k < d_linearity_tol ) {
312  if(d_rank == 0) {
313  std::cout << "linearly dependent sample! k = " << k << "\n";
314  std::cout << "d_linearity_tol = " << d_linearity_tol << "\n";
315  }
316  k = 0;
317  linearly_dependent_sample = true;
318  } else if ( d_num_samples >= d_max_basis_dimension || add_without_increase ) {
319  k = 0;
320  linearly_dependent_sample = true;
321  }
322  // Check to see if the "number of samples" (in IncrementalSVD and
323  // its subclasses, d_num_samples appears to be equal to the number
324  // of columns of the left singular vectors) is greater than or equal
325  // to the dimension of snapshot vectors. If so, then the vector
326  // addressed by the pointer u must be linearly dependent on the left
327  // singular vectors.
328  else if (d_num_samples >= d_total_dim) {
329  linearly_dependent_sample = true;
330  }
331  else {
332  linearly_dependent_sample = false;
333  }
334 
335  // Create Q.
336  double* Q;
337  constructQ(Q, l, k);
338  delete l;
339 
340  // Now get the singular value decomposition of Q.
341  Matrix* A;
342  Matrix* W;
343  Matrix* sigma;
344  bool result = svd(Q, A, sigma, W);
345 
346  // Done with Q.
347  delete [] Q;
348 
349  // If the svd was successful then add the sample. Otherwise clean up and
350  // return.
351  if (result) {
352 
353  // We need to add the sample if it is not linearly dependent or if it is
354  // linearly dependent and we are not skipping linearly dependent samples.
355  if ((linearly_dependent_sample && !d_skip_linearly_dependent) ) {
356  // This sample is linearly dependent and we are not skipping linearly
357  // dependent samples.
358  if(d_rank == 0) std::cout << "adding linearly dependent sample!\n";
359  addLinearlyDependentSample(A, W, sigma);
360  delete sigma;
361  }
362  else if (!linearly_dependent_sample) {
363  // This sample is not linearly dependent.
364 
365  // Compute j
366  Vector* j = u_vec.minus(basisl);
367  for (int i = 0; i < d_dim; ++i) {
368  j->item(i) /= k;
369  }
370 
371  // addNewSample will assign sigma to d_S hence it should not be
372  // deleted upon return.
373  addNewSample(j, A, W, sigma);
374  delete j;
375  }
376  delete basisl;
377  delete A;
378  delete W;
379 
380  // Compute the basis vectors.
381  computeBasis();
382  }
383  else {
384  delete basisl;
385  delete A;
386  delete W;
387  delete sigma;
388  }
389  return result;
390 }
391 
392 void
394  double*& Q,
395  const Vector* l,
396  double k)
397 {
398  CAROM_VERIFY(l != 0);
399  CAROM_VERIFY(l->dim() == numSamples());
400 
401  // Create Q.
402  Q = new double [(d_num_samples+1)*(d_num_samples+1)];
403 
404  // Fill Q in column major order.
405  int q_idx = 0;
406  for (int row = 0; row < d_num_samples; ++row) {
407  q_idx = row;
408  for (int col = 0; col < d_num_samples; ++col) {
409  if (row == col)
410  {
411  Q[q_idx] = d_S->item(col);
412  }
413  else
414  {
415  Q[q_idx] = 0.0;
416  }
417  q_idx += d_num_samples+1;
418  }
419  Q[q_idx] = l->item(row);
420  }
421  q_idx = d_num_samples;
422  for (int col = 0; col < d_num_samples; ++col) {
423  Q[q_idx] = 0.0;
424  q_idx += d_num_samples+1;
425  }
426  Q[q_idx] = k;
427 }
428 
429 bool
431  double* A,
432  Matrix*& U,
433  Matrix*& S,
434  Matrix*& V)
435 {
436  CAROM_VERIFY(A != 0);
437 
438  // Construct U, S, and V.
439  U = new Matrix(d_num_samples+1, d_num_samples+1, false);
440  S = new Matrix(d_num_samples+1, d_num_samples+1, false);
441  V = new Matrix(d_num_samples+1, d_num_samples+1, false);
442  for (int row = 0; row < d_num_samples+1; ++row) {
443  for (int col = 0; col < d_num_samples+1; ++col) {
444  S->item(row, col) = 0.0;
445  }
446  }
447 
448  // Use lapack's dgesdd Fortran function to perform the svd. As this is
449  // Fortran A and all the computed matrices are in column major order.
450  double* sigma = new double [d_num_samples+1];
451  char jobz = 'A';
452  int m = d_num_samples+1;
453  int n = d_num_samples+1;
454  int lda = d_num_samples+1;
455  int ldu = d_num_samples+1;
456  int ldv = d_num_samples+1;
457  int lwork = m*(4*m + 7);
458  double* work = new double [lwork];
459  int iwork[8*m];
460  int info;
461  dgesdd(&jobz,
462  &m,
463  &n,
464  A,
465  &lda,
466  sigma,
467  &U->item(0, 0),
468  &ldu,
469  &V->item(0, 0),
470  &ldv,
471  work,
472  &lwork,
473  iwork,
474  &info);
475  delete [] work;
476 
477  // If the svd succeeded, fill U and S. Otherwise clean up and return.
478  if (info == 0) {
479  // Place sigma into S.
480  for (int i = 0; i < d_num_samples+1; ++i) {
481  S->item(i, i) = sigma[i];
482  }
483  delete [] sigma;
484 
485  // U is column major order so convert it to row major order.
486  for (int row = 0; row < d_num_samples+1; ++row) {
487  for (int col = row+1; col < d_num_samples+1; ++col) {
488  double tmp = U->item(row, col);
489  U->item(row, col) = U->item(col, row);
490  U->item(col, row) = tmp;
491  }
492  }
493  /* if(d_update_right_SV) {
494  // V is column major order so convert it to row major order.
495  for (int row = 0; row < d_num_samples+1; ++row) {
496  for (int col = row+1; col < d_num_samples+1; ++col) {
497  double tmp = V->item(row, col);
498  V->item(row, col) = V->item(col, row);
499  V->item(col, row) = tmp;
500  }
501  }
502  }*/
503  }
504  else {
505  delete [] sigma;
506  }
507  return info == 0;
508 }
509 
510 double
512  const Matrix* m)
513 {
514  CAROM_ASSERT(m != 0);
515 
516  double result = 0.0;
517  if (d_num_samples > 1) {
518  int last_col = d_num_samples-1;
519  double tmp = 0.0;
520  int num_rows = m->numRows();
521  for (int i = 0; i < num_rows; ++i) {
522  tmp += m->item(i, 0) * m->item(i, last_col);
523  }
524  if (m->distributed() && d_size > 1) {
525  MPI_Allreduce(&tmp, &result, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
526  }
527  else {
528  result = tmp;
529  }
530  }
531  return result;
532 }
533 
534 }
void getInteger(const std::string &key, int &data)
Reads an integer associated with the supplied key from the currently open database file.
Definition: Database.h:181
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.
Definition: Database.h:94
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.
Definition: Database.cpp:38
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.
Definition: Matrix.h:177
int numColumns() const
Returns the number of columns in the Matrix. This method will return the same value from each process...
Definition: Matrix.h:220
Matrix * transposeMult(const Matrix &other) const
Multiplies the transpose of this Matrix with other and returns the product, reference version.
Definition: Matrix.h:642
Matrix * mult(const Matrix &other) const
Multiplies this Matrix with other and returns the product, reference version.
Definition: Matrix.h:283
const double & item(int row, int col) const
Const Matrix member access. Matrix data is stored in row-major format.
Definition: Matrix.h:1007
int numRows() const
Returns the number of rows of the Matrix on this processor.
Definition: Matrix.h:194
bool restore_state
If true the state of the incremental SVD will be restored when the object is created.
Definition: Options.h:323
double linearity_tol
The tolerance of the incremental SVD algorithm to determine whether or not a sample is linearly depen...
Definition: Options.h:271
bool save_state
If true the state of the incremental SVD will be written to disk when the object is deleted....
Definition: Options.h:317
int max_basis_dimension
The maximum dimension of the basis.
Definition: Options.h:235
Definition: SVD.h:29
Matrix * d_basis
The globalized basis vectors for the current time interval.
Definition: SVD.h:169
Matrix * d_snapshots
The globalized snapshot vectors for the current time interval.
Definition: SVD.h:209
Matrix * d_basis_right
The globalized right basis vectors for the current time interval.
Definition: SVD.h:177
bool isFirstSample() const
Returns true if the next sample will result in a new time interval.
Definition: SVD.h:138
Vector * d_S
The vector S which is small.
Definition: SVD.h:201
int d_num_samples
Number of samples stored for the current time interval.
Definition: SVD.h:151
Matrix * d_U
The matrix U which is large.
Definition: SVD.h:185
bool d_debug_algorithm
Flag to indicate if results of algorithm should be printed for debugging purposes.
Definition: SVD.h:215
const int d_dim
Dimension of the system.
Definition: SVD.h:146
Matrix * d_W
The matrix U which is large.
Definition: SVD.h:193
double norm() const
Form the norm of this.
Definition: Vector.cpp:242
Vector * minus(const Vector &other) const
Subtracts other and this and returns the result, reference version.
Definition: Vector.h:538
const double & item(int i) const
Const Vector member access.
Definition: Vector.h:656
double inner_product(const Vector &other) const
Inner product, reference form.
Definition: Vector.cpp:222
int dim() const
Returns the dimension of the Vector on this processor.
Definition: Vector.h:266