11 #include "linalg/Matrix.h"
12 #include "Utilities.h"
18 #include "STSampling.h"
30 void SampleTemporalIndices(
const Matrix* s_basis,
31 const Matrix* t_basis,
32 const int num_f_basis_vectors_used,
36 const int num_samples_req,
37 const bool excludeFinalTime)
39 CAROM_VERIFY(t_basis->distributed());
42 CAROM_VERIFY(0 < num_f_basis_vectors_used
43 && num_f_basis_vectors_used <= s_basis->numColumns()
44 && num_f_basis_vectors_used <= t_basis->numColumns());
45 CAROM_VERIFY(num_samples_req > 0);
46 const int num_basis_vectors =
47 std::min(num_f_basis_vectors_used, s_basis->numColumns());
48 const int num_samples = num_samples_req;
49 const int numExcluded = excludeFinalTime ? 1 : 0;
50 CAROM_VERIFY(num_samples <= t_basis->numRows() - numExcluded);
51 const int s_size = s_basis->numRows();
52 const int t_size = t_basis->numRows();
54 const int ns_mod_nr = num_samples % num_basis_vectors;
59 Matrix M(num_basis_vectors, num_basis_vectors,
62 std::set<int> samples;
64 std::vector<double> error(s_size * t_size);
65 std::vector<double> MZphi(num_basis_vectors);
68 for (
int s=0; s<s_size; ++s)
70 for (
int t=0; t<t_size; ++t)
71 error[t + (s*t_size)] = s_basis->item(s, 0) * t_basis->item(t, 0);
74 for (
int i=0; i<num_basis_vectors; ++i)
76 CAROM_VERIFY(samples.size() == ns);
83 M.setSize(ns * s_size, i);
86 for (
auto t : samples)
88 for (
int s = 0; s < s_size; ++s)
90 const int row = ti + (s*ns);
91 for (
int col = 0; col < i; ++col)
93 M.item(row, col) = s_basis->item(s, col) * t_basis->item(t, col);
101 M.transposePseudoinverse();
110 for (
int j = 0; j < i; ++j) MZphi[j] = 0.0;
113 for (
auto t : samples)
115 for (
int s = 0; s < s_size; ++s)
117 const int row = ti + (s*ns);
118 const double phi_i_s_t = s_basis->item(s, i) * t_basis->item(t,
121 for (
int col = 0; col < i; ++col)
123 MZphi[col] += M.item(row, col) * phi_i_s_t;
131 for (
int s=0; s<s_size; ++s)
133 for (
int t=0; t<t_size; ++t)
134 error[t + (s*t_size)] = s_basis->item(s, i) * t_basis->item(t, i);
138 for (
int s=0; s<s_size; ++s)
140 for (
int t=0; t<t_size; ++t)
142 for (
int j=0; j<i; ++j)
143 error[t + (s*t_size)] -= MZphi[j] * s_basis->item(s, j) * t_basis->item(t, j);
150 const int nsi = i < ns_mod_nr ? (num_samples / num_basis_vectors) + 1 :
151 num_samples / num_basis_vectors;
153 for (
int j=0; j<nsi; ++j)
157 double maxNorm = -1.0;
159 for (
int t=0; t<t_size - numExcluded; ++t)
161 std::set<int>::const_iterator found = samples.find(t);
162 if (found != samples.end())
166 for (
int s=0; s<s_size; ++s)
168 norm2 += error[t + (s*t_size)] * error[t + (s*t_size)];
172 double globalNorm2 = 0.0;
173 MPI_Allreduce(&norm2, &globalNorm2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
174 if (globalNorm2 > maxNorm)
176 maxNorm = globalNorm2;
181 samples.insert(tmax);
187 CAROM_VERIFY(num_samples == ns);
190 for (
auto t : samples)
203 void SampleSpatialIndices(
const Matrix* s_basis,
204 const Matrix* t_basis,
205 const int num_f_basis_vectors_used,
206 const int num_t_samples,
209 int* s_sampled_rows_per_proc,
210 Matrix& f_basis_sampled,
213 const int num_samples_req)
220 MPI_Datatype MaxRowType, oldtypes[2];
222 MPI_Aint offsets[2], extent, lower_bound;
225 oldtypes[0] = MPI_DOUBLE;
233 #if (MPI_VERSION == 1)
234 MPI_Type_extent(MPI_DOUBLE, &extent);
236 MPI_Type_get_extent(MPI_DOUBLE, &lower_bound, &extent);
240 oldtypes[1] = MPI_INT;
248 #if (MPI_VERSION == 1)
249 MPI_Type_struct(2, blockcounts, offsets, oldtypes, &MaxRowType);
251 MPI_Type_create_struct(2, blockcounts, offsets, oldtypes, &MaxRowType);
254 MPI_Type_commit(&MaxRowType);
258 MPI_Op_create((MPI_User_function*)RowInfoMax,
true, &RowInfoOp);
263 CAROM_VERIFY(0 < num_f_basis_vectors_used
264 && num_f_basis_vectors_used <= s_basis->numColumns()
265 && num_f_basis_vectors_used <= t_basis->numColumns());
266 CAROM_VERIFY(num_samples_req >
268 const int num_basis_vectors =
269 std::min(num_f_basis_vectors_used, s_basis->numColumns());
271 const int num_samples = num_samples_req;
272 CAROM_VERIFY(num_basis_vectors <= num_samples
273 && num_samples <= s_basis->numRows());
274 CAROM_VERIFY(num_samples == f_basis_sampled.numRows()
275 && num_basis_vectors == f_basis_sampled.numColumns());
276 CAROM_VERIFY(!f_basis_sampled.distributed());
277 const int s_size = s_basis->numRows();
278 const int t_size = t_basis->numRows();
280 const int ns_mod_nr = num_samples % num_basis_vectors;
285 Matrix M(num_basis_vectors, num_basis_vectors,
289 std::vector<std::set<int> > proc_sampled_f_row(num_procs);
290 std::vector<std::map<int, int> > proc_f_row_to_tmp_fs_row(num_procs);
291 int num_f_basis_cols =
292 f_basis_sampled.numColumns();
293 Matrix tmp_fs(f_basis_sampled.numRows(),
295 f_basis_sampled.distributed());
298 RowInfo f_bv_max_local, f_bv_max_global;
300 std::vector<double> error(s_size * t_size);
301 std::vector<double> MZphi(num_basis_vectors);
302 std::vector<double> sampled_row(num_basis_vectors);
305 for (
int s=0; s<s_size; ++s)
307 for (
int t=0; t<t_size; ++t)
308 error[t + (s*t_size)] = s_basis->item(s, 0) * t_basis->item(t, 0);
311 for (
int i=0; i<num_basis_vectors; ++i)
341 M.setSize(ns * num_t_samples, i);
342 for (
int k = 0; k < i; ++k)
343 for (
int j = 0; j < ns; ++j)
345 for (
int ti = 0; ti < num_t_samples; ++ti)
347 const int t = t_samples[ti];
348 const int row = ti + (j*num_t_samples);
350 M.item(row, k) = tmp_fs.item(j, k) * t_basis->item(t, k);
355 M.transposePseudoinverse();
359 for (
int j = 0; j < i; ++j) MZphi[j] = 0.0;
387 for (
int j = 0; j < ns; ++j)
389 for (
int ti = 0; ti < num_t_samples; ++ti)
391 const int t = t_samples[ti];
392 const int row = ti + (j*num_t_samples);
394 const double phi_i_s_t = s_basis->item(j, i) * t_basis->item(t,
397 for (
int k = 0; k < i; ++k)
398 MZphi[k] += M.item(row, k) * phi_i_s_t;
403 for (
int s=0; s<s_size; ++s)
405 for (
int t=0; t<t_size; ++t)
406 error[t + (s*t_size)] = s_basis->item(s, i) * t_basis->item(t, i);
410 for (
int s=0; s<s_size; ++s)
412 for (
int t=0; t<t_size; ++t)
414 for (
int j=0; j<i; ++j)
415 error[t + (s*t_size)] -= MZphi[j] * s_basis->item(s, j) * t_basis->item(t, j);
422 const int nsi = i < ns_mod_nr ? (num_samples / num_basis_vectors) + 1 :
423 num_samples / num_basis_vectors;
425 for (
int j=0; j<nsi; ++j)
430 double maxNorm = -1.0;
432 for (
int s=0; s<s_size; ++s)
434 std::set<int>::const_iterator found = proc_sampled_f_row[myid].find(s);
435 if (found != proc_sampled_f_row[myid].end())
439 for (
int t=0; t<t_size; ++t)
441 norm2 += error[t + (s*t_size)] * error[t + (s*t_size)];
454 f_bv_max_local.row_val = maxNorm;
455 f_bv_max_local.row = smax;
456 f_bv_max_local.proc = myid;
458 MPI_Allreduce(&f_bv_max_local, &f_bv_max_global, 1,
459 MaxRowType, RowInfoOp, MPI_COMM_WORLD);
461 proc_sampled_f_row[f_bv_max_global.proc].insert(f_bv_max_global.row);
462 proc_f_row_to_tmp_fs_row[f_bv_max_global.proc][f_bv_max_global.row] = ns + j;
465 if (f_bv_max_global.proc == myid) {
466 for (
int k = 0; k < num_basis_vectors; ++k) {
467 sampled_row[k] = s_basis->item(f_bv_max_global.row, k);
470 MPI_Bcast(sampled_row.data(), num_basis_vectors, MPI_DOUBLE,
471 f_bv_max_global.proc, MPI_COMM_WORLD);
473 for (
int k = 0; k < num_basis_vectors; ++k) {
474 tmp_fs.item(ns+j, k) = sampled_row[k];
481 CAROM_VERIFY(num_samples == ns);
486 for (
int i = 0; i < num_procs; ++i) {
487 std::set<int>& this_proc_sampled_f_row = proc_sampled_f_row[i];
488 std::map<int, int>& this_proc_f_row_to_tmp_fs_row =
489 proc_f_row_to_tmp_fs_row[i];
490 s_sampled_rows_per_proc[i] = this_proc_sampled_f_row.size();
491 for (std::set<int>::iterator j = this_proc_sampled_f_row.begin();
492 j != this_proc_sampled_f_row.end(); ++j) {
494 s_sampled_row[idx] = this_f_row;
495 int tmp_fs_row = this_proc_f_row_to_tmp_fs_row[this_f_row];
496 for (
int col = 0; col < num_f_basis_cols; ++col) {
497 f_basis_sampled.item(idx, col) = tmp_fs.item(tmp_fs_row, col);
503 CAROM_VERIFY(num_samples == idx);
506 MPI_Type_free(&MaxRowType);
507 MPI_Op_free(&RowInfoOp);
512 const int num_f_basis_vectors_used,
513 std::vector<int>& t_samples,
515 int* f_sampled_rows_per_proc,
519 const int num_t_samples_req,
520 const int num_s_samples_req,
521 const bool excludeFinalTime)
528 CAROM_VERIFY(s_basis->
numColumns() == num_f_basis_vectors_used
529 && t_basis->
numColumns() == num_f_basis_vectors_used);
532 CAROM_VERIFY(t_samples.size() == num_t_samples_req);
534 SampleTemporalIndices(s_basis, t_basis, num_f_basis_vectors_used,
536 myid, num_procs, num_t_samples_req, excludeFinalTime);
540 CAROM_VERIFY(s_basis_sampled.
numRows() == num_s_samples_req
541 && s_basis_sampled.
numColumns() == num_f_basis_vectors_used);
545 SampleSpatialIndices(s_basis, t_basis, num_f_basis_vectors_used,
548 t_samples.data(), f_sampled_row, f_sampled_rows_per_proc,
549 s_basis_sampled, myid, num_procs, num_s_samples_req);
552 void GetSampledSpaceTimeBasis(std::vector<int>
const& t_samples,
553 const Matrix* t_basis,
554 Matrix
const& s_basis_sampled,
555 Matrix& f_basis_sampled_inv)
557 const int num_s_samples = s_basis_sampled.numRows();
558 const int num_t_samples = t_samples.size();
562 CAROM_VERIFY(f_basis_sampled_inv.numRows() == num_t_samples * num_s_samples);
564 for (
int si=0; si<num_s_samples; ++si)
566 for (
int ti=0; ti<num_t_samples; ++ti)
568 const int row = ti + (si*num_t_samples);
569 const int t = t_samples[ti];
570 for (
int j=0; j<f_basis_sampled_inv.numColumns(); ++j)
571 f_basis_sampled_inv.item(row, j) = s_basis_sampled.item(si,
572 j) * t_basis->item(t, j);
577 f_basis_sampled_inv.transposePseudoinverse();
int numColumns() const
Returns the number of columns in the Matrix. This method will return the same value from each process...
int numRows() const
Returns the number of rows of the Matrix on this processor.
void SpaceTimeSampling(const Matrix *s_basis, const Matrix *t_basis, const int num_f_basis_vectors_used, std::vector< int > &t_samples, int *f_sampled_row, int *f_sampled_rows_per_proc, Matrix &s_basis_sampled, const int myid, const int num_procs, const int num_t_samples_req, const int num_s_samples_req, const bool excludeFinalTime)