diff --git a/examples/dmd/dg_advection.cpp b/examples/dmd/dg_advection.cpp
index 9e1d8b040..6bb905ac0 100644
--- a/examples/dmd/dg_advection.cpp
+++ b/examples/dmd/dg_advection.cpp
@@ -710,9 +710,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(U->Size());
     true_solution_u = U->GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (myid == 0)
     {
         std::cout << "Predicting solution using DMD" << std::endl;
@@ -755,9 +753,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_U.predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());
diff --git a/examples/dmd/dg_euler.cpp b/examples/dmd/dg_euler.cpp
index ae4f06885..30b596193 100644
--- a/examples/dmd/dg_euler.cpp
+++ b/examples/dmd/dg_euler.cpp
@@ -574,9 +574,7 @@ int main(int argc, char *argv[])
     Vector true_solution_e(u_block.GetBlock(3).Size());
     true_solution_e = u_block.GetBlock(3).GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (mpi.WorldRank() == 0)
     {
         std::cout << "Predicting density, momentum, and energy using DMD" << std::endl;
@@ -639,12 +637,12 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_dens = dmd_dens->predict(t_final);
     result_x_mom = dmd_x_mom->predict(t_final);
     result_y_mom = dmd_y_mom->predict(t_final);
     result_e = dmd_e->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_dens(result_dens->getData(), result_dens->dim());
diff --git a/examples/dmd/heat_conduction.cpp b/examples/dmd/heat_conduction.cpp
index baea9b2c7..3b229a87d 100644
--- a/examples/dmd/heat_conduction.cpp
+++ b/examples/dmd/heat_conduction.cpp
@@ -498,9 +498,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(u.Size());
     true_solution_u = u.GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
+    // 14. Predict using DMD.
     if (myid == 0)
     {
         std::cout << "Predicting temperature using DMD" << std::endl;
@@ -540,9 +538,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_u->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());
diff --git a/examples/dmd/heat_conduction_dmdc.cpp b/examples/dmd/heat_conduction_dmdc.cpp
index c1fa4d56e..91347c89a 100644
--- a/examples/dmd/heat_conduction_dmdc.cpp
+++ b/examples/dmd/heat_conduction_dmdc.cpp
@@ -534,9 +534,7 @@ int main(int argc, char *argv[])
     Vector true_solution_u(u.Size());
     true_solution_u = u.GetData();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMDc.
+    // 14. Predict using DMDc.
     if (myid == 0)
     {
         std::cout << "Predicting temperature using DMDc" << std::endl;
@@ -576,9 +574,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_u = dmd_u->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 15. Calculate the relative error between the DMDc final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());
diff --git a/examples/dmd/nonlinear_elasticity.cpp b/examples/dmd/nonlinear_elasticity.cpp
index 2c093c818..fce0eecb1 100644
--- a/examples/dmd/nonlinear_elasticity.cpp
+++ b/examples/dmd/nonlinear_elasticity.cpp
@@ -515,6 +515,12 @@ int main(int argc, char *argv[])
 
         if (last_step || (ti % windowNumSamples) == 0)
         {
+            // Calculate DMD modes
+            if (myid == 0 && rdim != -1 && ef != -1)
+            {
+                std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl;
+            }
+
             if (rdim != -1)
             {
                 if (myid == 0)
@@ -608,10 +614,10 @@ int main(int argc, char *argv[])
         w_gf.Save(ee_ofs);
     }
 
-    // 13. Calculate the DMD modes.
-    if (myid == 0 && rdim != -1 && ef != -1)
+    // 13. Predict using DMD.
+    if (myid == 0)
     {
-        std::cout << "Both rdim and ef are set. ef will be ignored." << std::endl;
+        std::cout << "Predicting position and velocity using DMD" << std::endl;
     }
 
     Vector true_solution_x(x_gf.GetTrueVector().Size());
@@ -620,14 +626,6 @@ int main(int argc, char *argv[])
     Vector true_solution_v(v_gf.GetTrueVector().Size());
     true_solution_v = v_gf.GetTrueVector();
 
-    dmd_prediction_timer.Start();
-
-    // 14. Predict the state at t_final using DMD.
-    if (myid == 0)
-    {
-        std::cout << "Predicting position and velocity using DMD" << std::endl;
-    }
-
     curr_window = 0;
     CAROM::Vector* result_x = dmd_x[curr_window]->predict(ts[0]);
     CAROM::Vector* result_v = dmd_v[curr_window]->predict(ts[0]);
@@ -687,12 +685,12 @@ int main(int argc, char *argv[])
         }
     }
 
-    dmd_prediction_timer.Stop();
-
+    dmd_prediction_timer.Start();
     result_x = dmd_x[curr_window]->predict(t_final);
     result_v = dmd_v[curr_window]->predict(t_final);
+    dmd_prediction_timer.Stop();
 
-    // 15. Calculate the relative error between the DMD final solution and the true solution.
+    // 14. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_x(result_x->getData(), result_x->dim());
     Vector diff_x(true_solution_x.Size());
     subtract(dmd_solution_x, true_solution_x, diff_x);
@@ -722,7 +720,7 @@ int main(int argc, char *argv[])
                dmd_prediction_timer.RealTime());
     }
 
-    // 16. Free the used memory.
+    // 15. Free the used memory.
     delete ode_solver;
     delete pmesh;
     delete result_x;
diff --git a/examples/dmd/wave_equation.cpp b/examples/dmd/wave_equation.cpp
index 52326e3e5..aa0f800cb 100644
--- a/examples/dmd/wave_equation.cpp
+++ b/examples/dmd/wave_equation.cpp
@@ -525,8 +525,7 @@ int main(int argc, char *argv[])
         dudt_gf.Save(osol);
     }
 
-    // 10. Predict the state at t_final using DMD.
-    dmd_prediction_timer.Start();
+    // 10. Predict using DMD.
     cout << "Predicting temperature using DMD" << endl;
     CAROM::Vector* result_u = nullptr;
     VisItDataCollection dmd_visit_dc("DMD_Wave_Equation", mesh);
@@ -564,8 +563,9 @@ int main(int argc, char *argv[])
             }
         }
     }
-    dmd_prediction_timer.Stop();
+    dmd_prediction_timer.Start();
     result_u = dmd_u[curr_window]->predict(t_final);
+    dmd_prediction_timer.Stop();
 
     // 11. Calculate the relative error between the DMD final solution and the true solution.
     Vector dmd_solution_u(result_u->getData(), result_u->dim());
diff --git a/examples/misc/combine_samples.cpp b/examples/misc/combine_samples.cpp
index 77987ab37..8ab5f19e1 100644
--- a/examples/misc/combine_samples.cpp
+++ b/examples/misc/combine_samples.cpp
@@ -53,15 +53,19 @@ int main(int argc, char* argv[])
     std::vector<std::string> sample_names;
     int snaps = 0;
     int dim   = 0;
+    int col_min = 1;
+    int col_max = 1e9;
     bool subtract_mean = false;
     bool subtract_offset = false;
+    std::string generator_filename = "total";
     std::string kind = "snapshot";
     std::string offset_file;
 
     bool offset_arg = false;
 
     if (argc >= 2) {
-        for (int i = 1; i < argc; i++) {
+        int i = 1;
+        while (i < argc) {
             if (!strcmp(argv[i], "basis") || !strcmp(argv[i], "-b")) {
                 if (rank==0) std::cout << "Argument " << i << " identified as basis or -b" <<
                                            std::endl;
@@ -80,9 +84,25 @@ int main(int argc, char* argv[])
                 if (rank==0) std::cout << "Will subtract mean" << std::endl;
                 subtract_mean = true;
             }
+            else if (!strcmp(argv[i], "file") || !strcmp(argv[i], "-f")) {
+                i += 1;
+                generator_filename = argv[i];
+                if (rank==0) std::cout << "Output prefix = " << generator_filename << std::endl;
+            }
+            else if (!strcmp(argv[i], "col_min") || !strcmp(argv[i], "-cmin")) {
+                i += 1;
+                col_min = std::stoi(argv[i]);
+                if (rank==0) std::cout << "First column to read = " << col_min << std::endl;
+            }
+            else if (!strcmp(argv[i], "col_max") || !strcmp(argv[i], "-cmax")) {
+                i += 1;
+                col_max = std::stoi(argv[i]);
+                if (rank==0) std::cout << "Last column to read = " << col_max << std::endl;
+            }
             else {
                 sample_names.push_back(argv[i]);
             }
+            i += 1;
         }
     }
     else {
@@ -112,7 +132,6 @@ int main(int argc, char* argv[])
     CAROM_VERIFY((snaps > 0) && (dim > 0));
 
     /*-- Load data from input files --*/
-    std::string generator_filename = "total";
     std::unique_ptr<CAROM::BasisGenerator> static_basis_generator;
     static_basis_generator.reset(new CAROM::BasisGenerator(
                                      CAROM::Options(dim, snaps).setMaxBasisDimension(snaps), false,
@@ -120,19 +139,22 @@ int main(int argc, char* argv[])
 
     if (rank==0) std::cout << "Loading data from " << kind << std::endl;
     for(const auto& sample_name: sample_names) {
-        static_basis_generator->loadSamples(sample_name, kind);
+        static_basis_generator->loadSampleRange(sample_name, kind, col_min-1,
+                                                col_max-1);
     }
 
     if (rank==0) std::cout << "Saving data uploaded as a snapshot matrix" <<
                                std::endl;
     static_basis_generator->writeSnapshot();
 
+    int rdim;
     if (!subtract_mean && !subtract_offset) {
         /*-- Compute SVD and save file --*/
         if (rank==0) std::cout << "Computing SVD" << std::endl;
         int rom_dim = static_basis_generator->getSpatialBasis()->numColumns();
         if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl;
         static_basis_generator->endSamples();
+        if (rank==0) static_basis_generator->finalSummary(1e-8, rdim);
     }
     else {
         /*-- load data from hdf5 file to find the mean and subtract it --*/
@@ -188,6 +210,7 @@ int main(int argc, char* argv[])
         int rom_dim = static_basis_generator2->getSpatialBasis()->numColumns();
         if (rank==0) std::cout << "U ROM Dimension: " << rom_dim << std::endl;
         static_basis_generator2->endSamples();
+        if (rank==0) static_basis_generator2->finalSummary(1e-8, rdim);
 
         static_basis_generator2 = nullptr;
     }
diff --git a/examples/prom/dg_advection_global_rom.cpp b/examples/prom/dg_advection_global_rom.cpp
index 4f28c57d9..2b62d8375 100644
--- a/examples/prom/dg_advection_global_rom.cpp
+++ b/examples/prom/dg_advection_global_rom.cpp
@@ -10,14 +10,21 @@
 
 //                       libROM MFEM Example: DG Advection (adapted from ex9p.cpp)
 //
+// =================================================================================
+//
 // Compile with: make dg_advection_global_rom
 //
-// For ROM (reproductive case):
+// Sample runs and results for reproductive case:
 //    dg_advection_global_rom -offline
 //    dg_advection_global_rom -merge -ns 1
 //    dg_advection_global_rom -online
 //
-// For ROM (global rom):
+// Output:
+//  Relative l2 error of ROM solution 0.0005641839849343575
+//
+// =================================================================================
+//
+// Sample runs and results for global ROM:
 // Offline phase: dg_advection_global_rom -offline -ff 1.0 -id 0
 //                dg_advection_global_rom -offline -ff 1.1 -id 1
 //                dg_advection_global_rom -offline -ff 1.2 -id 2
@@ -28,19 +35,10 @@
 //
 // Online phase:  dg_advection_global_rom -online -ff 1.15
 //
-// Sample runs:
-//    mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.005
-//    mpirun -np 4 dg_advection_global_rom -p 0 -dt 0.01
-//    mpirun -np 4 dg_advection_global_rom -p 1 -dt 0.005 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.002 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.02 -s 13 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.004 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 1 -rp 1 -dt 0.005 -tf 9
-//    mpirun -np 4 dg_advection_global_rom -p 3 -rp 2 -dt 0.0025 -tf 9 -vs 20
-//    mpirun -np 4 dg_advection_global_rom -p 0 -o 2 -rp 1 -dt 0.01 -tf 8
-//    mpirun -np 4 dg_advection_global_rom -p 0 -rs 2 -dt 0.005 -tf 2
-//    mpirun -np 4 dg_advection_global_rom -p 0 -rs 1 -o 2 -tf 2
-//    mpirun -np 3 dg_advection_global_rom -p 1 -rs 1 -rp 0 -dt 0.005 -tf 0.5
+// Output:
+//   Relative l2 error of ROM solution 0.0004333183604809453
+//
+// =================================================================================
 //
 // Description:  This example code solves the time-dependent advection equation
 //               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
@@ -128,10 +126,7 @@ class AIR_prec : public Solver
         delete AIR_solver;
         AIR_solver = new HypreBoomerAMG(A_s);
         AIR_solver->SetAdvectiveOptions(1, "", "FA");
-        AIR_solver->SetPrintLevel(
-            0);    // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
-        //    this mesh further in parallel to increase the resolution. Once the
-        //    parallel mesh is defined, the serial mesh can be deleted.
+        AIR_solver->SetPrintLevel(0);
         AIR_solver->SetMaxLevels(50);
     }
 
@@ -464,9 +459,6 @@ int main(int argc, char *argv[])
         {
             cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
         }
-        delete mesh;    // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
-        //    this mesh further in parallel to increase the resolution. Once the
-        //    parallel mesh is defined, the serial mesh can be deleted.
         return 3;
     }
 
diff --git a/examples/prom/mixed_nonlinear_diffusion.cpp b/examples/prom/mixed_nonlinear_diffusion.cpp
index 19ab15dd9..5d1e15f59 100644
--- a/examples/prom/mixed_nonlinear_diffusion.cpp
+++ b/examples/prom/mixed_nonlinear_diffusion.cpp
@@ -427,54 +427,6 @@ CAROM::Matrix* GetFirstColumns(const int N, const CAROM::Matrix* A)
 }
 
 // TODO: move this to the library?
-void BasisGeneratorFinalSummary(CAROM::BasisGenerator* bg,
-                                const double energyFraction, int & cutoff, const std::string cutoffOutputPath)
-{
-    const int rom_dim = bg->getSpatialBasis()->numColumns();
-    const CAROM::Vector* sing_vals = bg->getSingularValues();
-
-    MFEM_VERIFY(rom_dim <= sing_vals->dim(), "");
-
-    double sum = 0.0;
-    for (int sv = 0; sv < sing_vals->dim(); ++sv) {
-        sum += (*sing_vals)(sv);
-    }
-
-    vector<double> energy_fractions = {0.9999, 0.999, 0.99, 0.9};
-    bool reached_cutoff = false;
-
-    ofstream outfile(cutoffOutputPath);
-
-    double partialSum = 0.0;
-    for (int sv = 0; sv < sing_vals->dim(); ++sv) {
-        partialSum += (*sing_vals)(sv);
-        for (int i = energy_fractions.size() - 1; i >= 0; i--)
-        {
-            if (partialSum / sum > energy_fractions[i])
-            {
-                outfile << "For energy fraction: " << energy_fractions[i] << ", take first "
-                        << sv+1 << " of " << sing_vals->dim() << " basis vectors" << endl;
-                energy_fractions.pop_back();
-            }
-            else
-            {
-                break;
-            }
-        }
-
-        if (!reached_cutoff && partialSum / sum > energyFraction)
-        {
-            cutoff = sv+1;
-            reached_cutoff = true;
-        }
-    }
-
-    if (!reached_cutoff) cutoff = sing_vals->dim();
-    outfile << "Take first " << cutoff << " of " << sing_vals->dim() <<
-            " basis vectors" << endl;
-    outfile.close();
-}
-
 void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
                 std::string name)
 {
@@ -496,9 +448,10 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
     generator.endSamples(); // save the merged basis file
 
     int cutoff = 0;
-    BasisGeneratorFinalSummary(&generator, 0.9999, cutoff, "mergedSV_" + name);
+    generator.finalSummary(1e-4, cutoff, "mergedSV_" + name);
 }
 
+// TODO: move this to the library?
 const CAROM::Matrix* GetSnapshotMatrix(const int dimFOM, const int nparam,
                                        const int max_num_snapshots, std::string name)
 {
diff --git a/examples/prom/nonlinear_elasticity_global_rom.cpp b/examples/prom/nonlinear_elasticity_global_rom.cpp
index 7936d143c..911596035 100644
--- a/examples/prom/nonlinear_elasticity_global_rom.cpp
+++ b/examples/prom/nonlinear_elasticity_global_rom.cpp
@@ -307,68 +307,6 @@ CAROM::Matrix *GetFirstColumns(const int N, const CAROM::Matrix *A)
 }
 
 // TODO: move this to the library?
-void BasisGeneratorFinalSummary(CAROM::BasisGenerator *bg,
-                                const double energyFraction, int &cutoff,
-                                const std::string cutoffOutputPath)
-{
-    const int rom_dim = bg->getSpatialBasis()->numColumns();
-    const CAROM::Vector *sing_vals = bg->getSingularValues();
-
-    MFEM_VERIFY(rom_dim <= sing_vals->dim(), "");
-
-    double sum = 0.0;
-    for (int sv = 0; sv < sing_vals->dim(); ++sv)
-    {
-        sum += (*sing_vals)(sv);
-    }
-
-    vector<double> energy_fractions = {0.9999999, 0.999999, 0.99999, 0.9999, 0.999, 0.99, 0.9};
-    bool reached_cutoff = false;
-
-    ofstream outfile(cutoffOutputPath);
-
-    double partialSum = 0.0;
-    stringstream prec;
-    int ctr = 1;
-    char buffer[100];
-    for (int sv = 0; sv < sing_vals->dim(); ++sv)
-    {
-        partialSum += (*sing_vals)(sv);
-        for (int i = energy_fractions.size() - 1; i >= 0; i--)
-        {
-            if (partialSum / sum > energy_fractions[i])
-            {
-                // Format string
-                prec.str(std::string());
-                prec << "%." << ctr << "f";
-                sprintf(buffer, prec.str().c_str(), energy_fractions[i]);
-                ctr++;
-
-                // Output string
-                outfile << "For energy fraction: " << string(buffer) << ", take first "
-                        << sv + 1 << " of " << sing_vals->dim() << " basis vectors" << endl;
-                energy_fractions.pop_back();
-            }
-            else
-            {
-                break;
-            }
-        }
-
-        if (!reached_cutoff && partialSum / sum > energyFraction)
-        {
-            cutoff = sv + 1;
-            reached_cutoff = true;
-        }
-    }
-
-    if (!reached_cutoff)
-        cutoff = sing_vals->dim();
-    outfile << "Take first " << cutoff << " of " << sing_vals->dim() <<
-            " basis vectors" << endl;
-    outfile.close();
-}
-
 void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
                 std::string name)
 {
@@ -390,8 +328,7 @@ void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
     generator.endSamples(); // save the merged basis file
 
     int cutoff = 0;
-    BasisGeneratorFinalSummary(&generator, 0.9999999, cutoff,
-                               "mergedSV_" + name + ".txt");
+    generator.finalSummary(1e-7, cutoff, "mergedSV_" + name + ".txt");
 }
 
 const CAROM::Matrix *GetSnapshotMatrix(const int dimFOM, const int nparam,
diff --git a/lib/linalg/BasisGenerator.cpp b/lib/linalg/BasisGenerator.cpp
index 7c1925a78..c1a94a11b 100644
--- a/lib/linalg/BasisGenerator.cpp
+++ b/lib/linalg/BasisGenerator.cpp
@@ -21,6 +21,9 @@
 #include "svd/IncrementalSVDFastUpdate.h"
 #include "svd/IncrementalSVDBrand.h"
 
+#include <iomanip>
+#include <fstream>
+
 namespace CAROM {
 
 BasisGenerator::BasisGenerator(
@@ -161,10 +164,11 @@ BasisGenerator::takeSample(
 }
 
 void
-BasisGenerator::loadSamples(const std::string& base_file_name,
-                            const std::string& kind,
-                            int cut_off,
-                            Database::formats db_format)
+BasisGenerator::loadSampleRange(const std::string& base_file_name,
+                                const std::string& kind,
+                                int col_min,
+                                int col_max,
+                                Database::formats db_format)
 {
     CAROM_ASSERT(!base_file_name.empty());
     CAROM_VERIFY(kind == "basis" || kind == "snapshot");
@@ -186,10 +190,12 @@ BasisGenerator::loadSamples(const std::string& base_file_name,
 
     int num_rows = mat->numRows();
     int num_cols = mat->numColumns();
-    int max_cols = num_cols;
-    if (cut_off < num_cols) max_cols = cut_off;
+    if (col_min < 0) col_min = 0;
+    if (col_max > num_cols-1) col_max = num_cols-1;
+
+    CAROM_VERIFY(col_max >= col_min);
 
-    for (int j = 0; j < max_cols; j++) {
+    for (int j = col_min; j <= col_max; j++) {
         double* u_in = new double[num_rows];
         for (int i = 0; i < num_rows; i++) {
             if (kind == "basis") {
@@ -204,6 +210,15 @@ BasisGenerator::loadSamples(const std::string& base_file_name,
     }
 }
 
+void
+BasisGenerator::loadSamples(const std::string& base_file_name,
+                            const std::string& kind,
+                            int cutoff,
+                            Database::formats db_format)
+{
+    loadSampleRange(base_file_name, kind, 0, cutoff-1, db_format);
+}
+
 double
 BasisGenerator::computeNextSampleTime(
     double* u_in,
@@ -310,6 +325,79 @@ BasisGenerator::resetDt(
     }
 }
 
+void
+BasisGenerator::finalSummary(
+    const double energyFractionThreshold,
+    int & cutoff,
+    const std::string & cutoffOutputPath,
+    const int first_sv)
+{
+    const int rom_dim = getSpatialBasis()->numColumns();
+    const Vector* sing_vals = getSingularValues();
+
+    CAROM_VERIFY(rom_dim <= sing_vals->dim());
+
+    double sum = 0.0;
+    for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
+        sum += (*sing_vals)(sv);
+    }
+
+    int p = std::floor(-std::log10(energyFractionThreshold));
+    std::vector<double> energy_fractions(p);
+
+    for (int i = 0; i < p; ++i) {
+        energy_fractions[i] = 1 - std::pow(10, -1 - i);
+    }
+
+    cutoff = first_sv;
+    bool reached_cutoff = false;
+    double partialSum = 0.0;
+    int count = 0;
+
+    std::ostream* output_stream;
+
+    if (!cutoffOutputPath.empty()) {
+        output_stream = new std::ofstream(cutoffOutputPath);
+    } else {
+        output_stream = &std::cout;
+    }
+
+    for (int sv = first_sv; sv < sing_vals->dim(); ++sv) {
+        partialSum += (*sing_vals)(sv);
+        for (int i = count; i < p; ++i)
+        {
+            if (partialSum / sum > 1.0 - std::pow(10, -1 - i))
+            {
+                *output_stream << "For energy fraction: 0.";
+                for (int j = 0; j < i+1; ++j) *output_stream << "9";
+                *output_stream << ", take first " << sv+1 << " of "
+                               << sing_vals->dim() << " basis vectors" << std::endl;
+                count += 1;
+            }
+            else
+            {
+                break;
+            }
+        }
+        if (!reached_cutoff && partialSum / sum > 1.0 - energyFractionThreshold)
+        {
+            cutoff = sv+1;
+            reached_cutoff = true;
+        }
+    }
+
+    if (!reached_cutoff) cutoff = sing_vals->dim();
+    *output_stream << std::fixed << std::setprecision(p+1);
+    *output_stream << "For energy fraction: " << 1.0 - energyFractionThreshold <<
+                   ", take first "
+                   << cutoff << " of " << sing_vals->dim() << " basis vectors" << std::endl;
+
+    if (!cutoffOutputPath.empty()) {
+        static_cast<std::ofstream*>(output_stream)->close();
+        delete output_stream;
+    }
+}
+
 BasisGenerator::~BasisGenerator()
 {
     if (d_basis_writer) {
diff --git a/lib/linalg/BasisGenerator.h b/lib/linalg/BasisGenerator.h
index ebad59495..2857db8d3 100644
--- a/lib/linalg/BasisGenerator.h
+++ b/lib/linalg/BasisGenerator.h
@@ -142,20 +142,39 @@ class BasisGenerator
         }
     }
 
+    /**
+     * @brief Load previously saved sample (basis or state)
+     *        within a column range.
+     *
+     * @param[in] base_file_name The base part of the name of the files
+     *                           holding the basis/snapshot vectors.
+     * @param[in] kind A string equal to "basis" or "snapshot", representing
+     *                 which kind of data to load.
+     * @param[in] col_min The first basis/snapshot vector to read.
+     * @param[in] col_max The last basis/snapshot vector to read.
+     * @param[in] db_format Format of the file to read.
+     */
+    void
+    loadSampleRange(const std::string& base_file_name,
+                    const std::string& kind  = "basis",
+                    int col_min = 0,
+                    int col_max = 1e9,
+                    Database::formats db_format = Database::HDF5);
+
     /**
      * @brief Load previously saved sample (basis or state).
      *
      * @param[in] base_file_name The base part of the name of the files
-     *                           holding the basis / snapshot vectors.
+     *                           holding the basis/snapshot vectors.
      * @param[in] kind A string equal to "basis" or "snapshot", representing
      *                 which kind of data to load.
-     * @param[in] cut_off The maximum number of bases or snapshots to read.
+     * @param[in] cutoff The maximum number of basis/snapshot vectors to read.
      * @param[in] db_format Format of the file to read.
      */
     void
     loadSamples(const std::string& base_file_name,
                 const std::string& kind  = "basis",
-                int cut_off = 1e9,
+                int cutoff = 1e9,
                 Database::formats db_format = Database::HDF5);
 
     /**
@@ -263,6 +282,21 @@ class BasisGenerator
         return d_svd->getNumSamples();
     }
 
+    /**
+     * @brief Prints the summary of recommended numbers of basis vectors.
+     *
+     * @param[in] energyFractionThreshold   Energy Fraction threshold
+     *                                      (energy fraction = 1.0 - energyFractionThreshold).
+     * @param[in] cutoff                    Number of basis vectors selected.
+     * @param[in] cutoffOutputPath          Path of the summary file.
+     * @param[in] first_sv                  First singular vector in the calculaton of energy.
+     */
+    void finalSummary(
+        const double energyFractionThreshold,
+        int & cutoff,
+        const std::string & cutoffOutputPath = "",
+        const int first_sv = 0);
+
 protected:
     /**
      * @brief Writer of basis vectors.