libROM  v1.0
Data-driven physical simulation library
StaticSVD.cpp
1 
11 // Description: A class implementing interface of SVD for the static SVD
12 // algorithm.
13 
14 #include "StaticSVD.h"
15 
16 #include "mpi.h"
17 #include "linalg/scalapack_wrapper.h"
18 
19 #include <limits.h>
20 
21 #include <stdio.h>
22 #include <string.h>
23 
24 /* Use automatically detected Fortran name-mangling scheme */
25 #define dgesdd CAROM_FC_GLOBAL(dgesdd, DGESDD)
26 
27 extern "C" {
28  void dgesdd(char*, int*, int*, double*, int*,
29  double*, double*, int*, double*, int*,
30  double*, int*, int*, int*);
31 }
32 
33 namespace CAROM {
34 
35 StaticSVD::StaticSVD(
36  Options options) :
37  SVD(options),
38  d_samples(new SLPK_Matrix), d_factorizer(new SVDManager),
39  d_basis_is_current(false),
40  d_max_basis_dimension(options.max_basis_dimension),
41  d_singular_value_tol(options.singular_value_tol),
42  d_preserve_snapshot(options.static_svd_preserve_snapshot)
43 {
44  // Get the rank of this process, and the number of processors.
45  int mpi_init;
46  MPI_Initialized(&mpi_init);
47  if (mpi_init == 0) {
48  MPI_Init(nullptr, nullptr);
49  }
50 
51  MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
52  MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs);
53 
55  /* TODO: Try doing this more intelligently and see if it makes a difference */
57  d_npcol = 1;
59  if (d_total_dim % d_nprow != 0) {
60  d_blocksize += 1;
61  }
62 
63  initialize_matrix(d_samples.get(), d_total_dim, d_max_num_samples,
64  d_nprow, d_npcol, d_blocksize, d_blocksize); // TODO: should nb = 1?
65  d_factorizer->A = nullptr;
66 }
67 
69 {
72 }
73 
75 {
76  if (d_samples) {
77  free_matrix_data(d_samples.get());
78  release_context(d_samples.get());
79  }
80 }
81 
83 {
84  if (d_factorizer->A != nullptr) {
85  free(d_factorizer->S);
86  d_factorizer->S = nullptr;
87  if (d_factorizer->U != nullptr) {
88  free_matrix_data(d_factorizer->U);
89  }
90  free(d_factorizer->U);
91  d_factorizer->U = nullptr;
92  if (d_factorizer->V != nullptr) {
93  free_matrix_data(d_factorizer->V);
94  }
95  free(d_factorizer->V);
96  d_factorizer->V = nullptr;
97  }
98 }
99 
100 bool
102  double* u_in,
103  bool add_without_increase)
104 {
105  CAROM_VERIFY(u_in != 0);
106  CAROM_NULL_USE(add_without_increase);
107  CAROM_VERIFY(0 <= d_num_samples);
108  CAROM_VERIFY(d_num_samples < d_max_num_samples);
109 
110  // Check the u_in is not non-zero.
111  Vector u_vec(u_in, d_dim, true);
112  if (u_vec.norm() == 0.0) {
113  return false;
114  }
115 
116  if (isFirstSample()) {
118  d_num_samples = 0;
119  d_basis = nullptr;
120  d_basis_right = nullptr;
121  // Set the N in the global matrix so BLACS won't complain.
123  }
124  broadcast_sample(u_in);
125  ++d_num_samples;
126 
127  // Build snapshot matrix before SVD is computed
128  //d_snapshots = new Matrix(d_dim, d_num_samples, false);
129  //for (int rank = 0; rank < d_num_procs; ++rank) {
130  // int nrows = d_dims[static_cast<unsigned>(rank)];
131  // int firstrow = d_istarts[static_cast<unsigned>(rank)] + 1;
132  // gather_transposed_block(&d_snapshots->item(0, 0), d_samples.get(),
133  // firstrow, 1, nrows,
134  // d_num_samples, rank);
135  // }
136  d_basis_is_current = false;
137  return true;
138 }
139 
140 const Matrix*
142 {
143  // If this basis is for the last time interval then it may not be up to date
144  // so recompute it.
145  if (!isBasisCurrent()) {
146  delete d_basis;
147  d_basis = nullptr;
148  delete d_basis_right;
149  d_basis_right = nullptr;
150  delete d_U;
151  d_U = nullptr;
152  delete d_S;
153  d_S = nullptr;
154  delete d_W;
155  d_W = nullptr;
156  computeSVD();
157  }
158  else {
159  CAROM_ASSERT(d_basis != 0);
160  }
161  CAROM_ASSERT(isBasisCurrent());
162  return d_basis;
163 }
164 
165 const Matrix*
167 {
168  // If this basis is for the last time interval then it may not be up to date
169  // so recompute it.
170  if (!isBasisCurrent()) {
171  delete d_basis;
172  d_basis = nullptr;
173  delete d_basis_right;
174  d_basis_right = nullptr;
175  delete d_U;
176  d_U = nullptr;
177  delete d_S;
178  d_S = nullptr;
179  delete d_W;
180  d_W = nullptr;
181  computeSVD();
182  }
183  else {
184  CAROM_ASSERT(d_basis_right != 0);
185  }
186  CAROM_ASSERT(isBasisCurrent());
187  return d_basis_right;
188 }
189 
190 const Vector*
192 {
193  // If these singular values are for the last time interval then they may not
194  // be up to date so recompute them.
195  if (!isBasisCurrent()) {
196  delete d_basis;
197  d_basis = nullptr;
198  delete d_basis_right;
199  d_basis = nullptr;
200  delete d_U;
201  d_U = nullptr;
202  delete d_S;
203  d_S = nullptr;
204  delete d_W;
205  d_W = nullptr;
206  computeSVD();
207  }
208  else {
209  CAROM_ASSERT(d_S != 0);
210  }
211  CAROM_ASSERT(isBasisCurrent());
212  return d_S;
213 }
214 
215 const Matrix*
217 {
218  if ((!d_preserve_snapshot) && (isBasisCurrent()) && (d_basis != 0))
219  CAROM_ERROR("StaticSVD: snapshot matrix is modified after computeSVD."
220  " To preserve the snapshots, set Options::static_svd_preserve_snapshot to be true!\n");
221 
222  if (d_snapshots) delete d_snapshots;
223  d_snapshots = new Matrix(d_dim, d_num_samples, true);
224 
225  for (int rank = 0; rank < d_num_procs; ++rank) {
226  int nrows = d_dims[static_cast<unsigned>(rank)];
227  int firstrow = d_istarts[static_cast<unsigned>(rank)] + 1;
228  gather_transposed_block(&d_snapshots->item(0, 0), d_samples.get(),
229  firstrow, 1, nrows,
230  d_num_samples, rank);
231  }
232 
233  CAROM_ASSERT(d_snapshots != 0);
234  return d_snapshots;
235 }
236 
237 void
239 {
240  // pointer for snapshot matrix or its transpose.
241  SLPK_Matrix *snapshot_matrix = NULL;
242 
243  // Finalizes the column size of d_samples.
245 
246  // transpose matrix if sample size > dimension.
247  bool transpose = d_total_dim < d_num_samples;
248 
249  if (transpose) {
250  // create a transposed matrix if sample size > dimension.
251  snapshot_matrix = new SLPK_Matrix;
252  int d_blocksize_tr = d_num_samples / d_nprow;
253  if (d_num_samples % d_nprow != 0) {
254  d_blocksize_tr += 1;
255  }
256 
257  initialize_matrix(snapshot_matrix, d_num_samples, d_total_dim,
258  d_nprow, d_npcol, d_blocksize_tr, d_blocksize_tr);
259 
260  for (int rank = 0; rank < d_num_procs; ++rank) {
261  transpose_submatrix(snapshot_matrix, 1,
262  d_istarts[static_cast<unsigned>(rank)]+1,
263  d_samples.get(), d_istarts[static_cast<unsigned>(rank)]+1, 1,
264  d_dims[static_cast<unsigned>(rank)], d_num_samples);
265  }
266  }
267  else {
268  // use d_samples if sample size <= dimension.
269  // SVD operation does not preserve the original snapshot matrix.
270  if (d_preserve_snapshot)
271  {
272  snapshot_matrix = new SLPK_Matrix;
273  initialize_matrix(snapshot_matrix, d_total_dim, d_num_samples,
275  copy_matrix(snapshot_matrix, 1, 1, d_samples.get(), 1, 1, d_total_dim,
276  d_num_samples);
277  }
278  else
279  // use the original snapshot matrix, which will be modified.
280  snapshot_matrix = d_samples.get();
281  }
282 
283  // This block does the actual ScaLAPACK call to do the factorization.
285  svd_init(d_factorizer.get(), snapshot_matrix);
286  d_factorizer->dov = 1;
287  factorize(d_factorizer.get());
288 
289  // Compute how many basis vectors we will actually return.
290  int sigma_cutoff = 0, hard_cutoff = d_num_samples;
291  if (d_singular_value_tol == 0) {
292  sigma_cutoff = std::numeric_limits<int>::max();
293  } else {
294  for (int i = 0; i < d_num_samples; ++i) {
295  if (d_factorizer->S[i] / d_factorizer->S[0] > d_singular_value_tol) {
296  sigma_cutoff += 1;
297  } else {
298  break;
299  }
300  }
301  }
302  if (d_max_basis_dimension != -1 && d_max_basis_dimension < hard_cutoff) {
303  hard_cutoff = d_max_basis_dimension;
304  }
305  int ncolumns = hard_cutoff < sigma_cutoff ? hard_cutoff : sigma_cutoff;
306  if (transpose) ncolumns = (ncolumns > d_total_dim) ? d_total_dim : ncolumns;
307  CAROM_VERIFY(ncolumns >= 0);
308 
309  // Allocate the appropriate matrices and gather their elements.
310  d_S = new Vector(ncolumns, false);
311  {
312  CAROM_VERIFY(ncolumns >= 0);
313  unsigned nc = static_cast<unsigned>(ncolumns);
314  memset(&d_S->item(0), 0, nc*sizeof(double));
315  }
316 
317  d_basis = new Matrix(d_dim, ncolumns, true);
318  d_basis_right = new Matrix(d_num_samples, ncolumns, false);
319  for (int rank = 0; rank < d_num_procs; ++rank) {
320  if (transpose) {
321  // V is computed in the transposed order so no reordering necessary.
322  gather_block(&d_basis->item(0, 0), d_factorizer->V,
323  1, d_istarts[static_cast<unsigned>(rank)]+1,
324  ncolumns, d_dims[static_cast<unsigned>(rank)], rank);
325 
326  // gather_transposed_block does the same as gather_block, but transposes
327  // it; here, it is used to go from column-major to row-major order.
328  gather_transposed_block(&d_basis_right->item(0, 0), d_factorizer->U,
329  1, 1, d_num_samples, ncolumns, rank);
330 
331  }
332  else {
333  // gather_transposed_block does the same as gather_block, but transposes
334  // it; here, it is used to go from column-major to row-major order.
335  gather_transposed_block(&d_basis->item(0, 0), d_factorizer->U,
336  d_istarts[static_cast<unsigned>(rank)]+1,
337  1, d_dims[static_cast<unsigned>(rank)],
338  ncolumns, rank);
339  // V is computed in the transposed order so no reordering necessary.
340  gather_block(&d_basis_right->item(0, 0), d_factorizer->V, 1, 1,
341  ncolumns, d_num_samples, rank);
342  }
343  }
344  for (int i = 0; i < ncolumns; ++i)
345  d_S->item(i) = d_factorizer->S[static_cast<unsigned>(i)];
346  d_basis_is_current = true;
347 
348  if (d_debug_algorithm) {
349  if (d_rank == 0) {
350  printf("Distribution of sampler's A and U:\n");
351  }
352  print_debug_info(d_samples.get());
353  MPI_Barrier(MPI_COMM_WORLD);
354 
355  if (d_rank == 0) {
356  printf("Distribution of sampler's V:\n");
357  }
358  print_debug_info(d_factorizer->V);
359  MPI_Barrier(MPI_COMM_WORLD);
360 
361  if (d_rank == 0) {
362  printf("Computed singular values: ");
363  for (int i = 0; i < ncolumns; ++i)
364  printf("%8.4E ", d_factorizer->S[i]);
365  printf("\n");
366  }
367  }
368 
369  // Delete the snapshot matrix.
370  if ((transpose) || (d_preserve_snapshot))
371  {
372  free_matrix_data(snapshot_matrix);
373  delete snapshot_matrix;
374  }
375 }
376 
377 void
378 StaticSVD::broadcast_sample(const double* u_in)
379 {
380  for (int rank = 0; rank < d_num_procs; ++rank) {
381  scatter_block(d_samples.get(), d_istarts[static_cast<unsigned>(rank)]+1,
382  d_num_samples+1, u_in, d_dims[static_cast<unsigned>(rank)],
383  1, rank);
384  }
385 }
386 
387 void
389 {
390  d_dims.resize(static_cast<unsigned>(d_num_procs));
391  MPI_Allgather(&d_dim, 1, MPI_INT, d_dims.data(), 1, MPI_INT, MPI_COMM_WORLD);
392  d_total_dim = 0;
393  d_istarts = std::vector<int>(static_cast<unsigned>(d_num_procs), 0);
394 
395  for (unsigned i = 0; i < static_cast<unsigned>(d_num_procs); ++i) {
396  d_total_dim += d_dims[i];
397  if (i > 0) {
398  d_istarts[i] = d_istarts[i-1] + d_dims[i-1];
399  }
400  }
401 }
402 
403 }
const double & item(int row, int col) const
Const Matrix member access. Matrix data is stored in row-major format.
Definition: Matrix.h:1015
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
const int d_max_num_samples
The maximum number of samples.
Definition: SVD.h:161
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 d_singular_value_tol
The tolerance for singular values below which to drop vectors.
Definition: StaticSVD.h:205
virtual bool takeSample(double *u_in, bool add_without_increase=false)
Collect the new sample, u_in at the supplied time.
Definition: StaticSVD.cpp:101
int d_rank
The rank of the process this object belongs to.
Definition: StaticSVD.h:150
int d_nprow
The number of processor rows in the grid.
Definition: StaticSVD.h:179
virtual void computeSVD()
Gathers samples from all other processors to form complete sample of system and computes the SVD.
Definition: StaticSVD.cpp:238
int d_npcol
The number of processor columns in the grid.
Definition: StaticSVD.h:184
void delete_factorizer()
Delete the factorizer from ScaLAPACK.
Definition: StaticSVD.cpp:82
void get_global_info()
Get the system's total row dimension and where my rows sit in the matrix.
Definition: StaticSVD.cpp:388
void delete_samples()
Delete the samples from ScaLAPACK.
Definition: StaticSVD.cpp:74
void broadcast_sample(const double *u_in)
Broadcast the sample to all processors.
Definition: StaticSVD.cpp:378
bool d_basis_is_current
Flag to indicate if the basis vectors for the current time interval are up to date.
Definition: StaticSVD.h:145
virtual const Vector * getSingularValues()
Returns the singular values for the current time interval.
Definition: StaticSVD.cpp:191
std::vector< int > d_istarts
The starting row (0-based) of the matrix that each process owns. Stored to avoid an MPI operation to ...
Definition: StaticSVD.h:162
int d_total_dim
The total dimension of the system (row dimension)
Definition: StaticSVD.h:174
int d_max_basis_dimension
The max number of basis vectors to return.
Definition: StaticSVD.h:200
std::unique_ptr< SLPK_Matrix > d_samples
Current samples of the system.
Definition: StaticSVD.h:134
std::unique_ptr< SVDManager > d_factorizer
Factorization manager object used to compute the SVD.
Definition: StaticSVD.h:139
int d_blocksize
The block size used internally for computing the SVD.
Definition: StaticSVD.h:189
virtual const Matrix * getSpatialBasis()
Returns the basis vectors for the current time interval as a Matrix.
Definition: StaticSVD.cpp:141
virtual const Matrix * getSnapshotMatrix()
Returns the snapshot matrix for the current time interval.
Definition: StaticSVD.cpp:216
bool isBasisCurrent()
Checks if the basis vectors are computed from the current snapshot.
Definition: StaticSVD.h:126
virtual const Matrix * getTemporalBasis()
Returns the temporal basis vectors for the current time interval as a Matrix.
Definition: StaticSVD.cpp:166
int d_num_procs
The number of processors being run on.
Definition: StaticSVD.h:155
std::vector< int > d_dims
The number of rows that each process owns. Stored to avoid an MPI operation to get this operation eve...
Definition: StaticSVD.h:169
double norm() const
Form the norm of this.
Definition: Vector.cpp:242
const double & item(int i) const
Const Vector member access.
Definition: Vector.h:656