A fix for StaticSVD::ComputeSVD to preserve snapshots

lib/linalg/Options.h

@@ -345,6 +345,11 @@ class Options
      *        incremental SVD algorithm.
      */
     double max_sampling_time_step_scale = 5.0;
+
+    /**
+     * @brief Option to preserve snapshot in StaticSVD::computeSVD.
+     */
+    bool static_svd_preserve_snapshot = false;
 };
 
 }

lib/linalg/svd/StaticSVD.cpp

@@ -38,7 +38,8 @@ StaticSVD::StaticSVD(
     d_samples(new SLPK_Matrix), d_factorizer(new SVDManager),
     d_this_interval_basis_current(false),
     d_max_basis_dimension(options.max_basis_dimension),
-    d_singular_value_tol(options.singular_value_tol)
+    d_singular_value_tol(options.singular_value_tol),
+    d_preserve_snapshot(options.static_svd_preserve_snapshot)
 {
     // Get the rank of this process, and the number of processors.
     int mpi_init;
@@ -234,6 +235,9 @@ StaticSVD::getSingularValues()
 const Matrix*
 StaticSVD::getSnapshotMatrix()
 {
+    if ((!d_preserve_snapshot) && (thisIntervalBasisCurrent()) && (d_basis != 0))
+        CAROM_ERROR("StaticSVD: snapshot matrix is modified after computeSVD."
+                    " To preserve the snapshots, set Options::static_svd_preserve_snapshot to be true!\n");
 
     if (d_snapshots) delete d_snapshots;
     d_snapshots = new Matrix(d_dim, d_num_samples, false);
@@ -265,7 +269,6 @@ StaticSVD::computeSVD()
     if (transpose) {
         // create a transposed matrix if sample size > dimension.
         snapshot_matrix = new SLPK_Matrix;
-
         int d_blocksize_tr = d_num_samples / d_nprow;
         if (d_num_samples % d_nprow != 0) {
             d_blocksize_tr += 1;
@@ -283,7 +286,18 @@ StaticSVD::computeSVD()
     }
     else {
         // use d_samples if sample size <= dimension.
-        snapshot_matrix = d_samples.get();
+        // SVD operation does not preserve the original snapshot matrix.
+        if (d_preserve_snapshot)
+        {
+            snapshot_matrix = new SLPK_Matrix;
+            initialize_matrix(snapshot_matrix, d_total_dim, d_num_samples,
+                              d_nprow, d_npcol, d_blocksize, d_blocksize);
+            copy_matrix(snapshot_matrix, 1, 1, d_samples.get(), 1, 1, d_total_dim,
+                        d_num_samples);
+        }
+        else
+            // use the original snapshot matrix, which will be modified.
+            snapshot_matrix = d_samples.get();
     }
 
     // This block does the actual ScaLAPACK call to do the factorization.
@@ -372,8 +386,9 @@ StaticSVD::computeSVD()
         }
     }
 
-    if (transpose) {
-        // Delete the transposed snapshot matrix.
+    // Delete the snapshot matrix.
+    if ((transpose) || (d_preserve_snapshot))
+    {
         free_matrix_data(snapshot_matrix);
         delete snapshot_matrix;
     }

lib/linalg/svd/StaticSVD.h

@@ -244,6 +244,11 @@ class StaticSVD : public SVD
     StaticSVD&
     operator = (
         const StaticSVD& rhs);
+
+    /**
+     * @brief option to preserve snapshot in computeSVD.
+     */
+    bool d_preserve_snapshot;
 };
 
 }