diff --git a/CMakeLists.txt b/CMakeLists.txt
index b8a0a5a..3f4b5d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,6 +5,8 @@ project(_pylibROM)
 set(CMAKE_BUILD_TYPE Debug)
 set(PYBIND11_FINDPYTHON ON)
 
+option (USE_MFEM "Build pylibROM with MFEM" OFF)
+
 #=================== ScaLAPACK (optional) ==================
 option(BUILD_SCALAPACK "Build static ScaLAPACK for libROM" OFF)
 
@@ -50,14 +52,18 @@ if (BUILD_LIBROM)
 	# )
 	# add_custom_target(RUN_LIBROM_BUILD ALL DEPENDS LIBROM_BUILD)
 
-	ExternalProject_Add(
-        	libROM
-        	SOURCE_DIR ${LIBROM_SCRIPTS_DIR}
-        	CONFIGURE_COMMAND ""
-        	BINARY_DIR ${LIBROM_DIR}
-        	BUILD_COMMAND ${LIBROM_SCRIPTS_DIR}/compile.sh -m -g -t ${LIBROM_DIR}/cmake/toolchains/simple.cmake
-        	INSTALL_COMMAND ""
-	)	
+    set(LIBROM_BUILD_CMD "${LIBROM_SCRIPTS_DIR}/compile.sh -t ${LIBROM_DIR}/cmake/toolchains/simple.cmake" CACHE STRING  "Command used to build libROM and dependencies")
+    if (USE_MFEM)
+        set(LIBROM_BUILD_CMD "${LIBROM_BUILD_CMD} -m -g")
+    endif()
+	# ExternalProject_Add(
+    #     	libROM
+    #     	SOURCE_DIR ${LIBROM_SCRIPTS_DIR}
+    #     	CONFIGURE_COMMAND ""
+    #     	BINARY_DIR ${LIBROM_DIR}
+    #         BUILD_COMMAND ${LIBROM_BUILD_CMD}
+    #         INSTALL_COMMAND ""
+    # )
 	message("Building libROM dependency...")
 endif(BUILD_LIBROM)
 
@@ -72,32 +78,34 @@ execute_process(COMMAND python3 -c "import mpi4py; print(mpi4py.get_include())"
 # # TODO(kevin): We do not bind mfem-related functions until we figure out how to type-cast SWIG Object.
 # #              Until then, mfem-related functions need to be re-implemented on python-end, using PyMFEM.
 
-find_library(MFEM mfem
-	"$ENV{MFEM_DIR}/lib"
-	"$ENV{MFEM_DIR}"
-	"${LIBROM_DIR}/dependencies/mfem")
-find_library(HYPRE HYPRE
-	"$ENV{HYPRE_DIR}/lib"
-	"${LIBROM_DIR}/dependencies/hypre/src/hypre/lib")
-find_library(PARMETIS parmetis
-	"$ENV{PARMETIS_DIR}/lib"
-	"$ENV{PARMETIS_DIR}/build/lib/libparmetis"
-	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libparmetis")
-find_library(METIS metis
-	"$ENV{METIS_DIR}/lib"
-	"$ENV{PARMETIS_DIR}/build/lib/libmetis"
-	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libmetis")
-find_path(MFEM_INCLUDES mfem.hpp
-	"$ENV{MFEM_DIR}/include"
-	"$ENV{MFEM_DIR}"
-	"${LIBROM_DIR}/dependencies/mfem")
-find_path(HYPRE_INCLUDES HYPRE.h
-	"$ENV{HYPRE_DIR}/include"
-	"${LIBROM_DIR}/dependencies/hypre/src/hypre/include")
-find_path(PARMETIS_INCLUDES metis.h
-	"$ENV{PARMETIS_DIR}/metis/include"
-	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/metis/include")
-
+if (USE_MFEM)
+    find_library(MFEM mfem
+        "$ENV{MFEM_DIR}/lib"
+        "$ENV{MFEM_DIR}"
+        "${LIBROM_DIR}/dependencies/mfem")
+    find_library(HYPRE HYPRE
+        "$ENV{HYPRE_DIR}/lib"
+        "${LIBROM_DIR}/dependencies/hypre/src/hypre/lib")
+    find_library(PARMETIS parmetis
+        "$ENV{PARMETIS_DIR}/lib"
+        "$ENV{PARMETIS_DIR}/build/lib/libparmetis"
+        "${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libparmetis")
+    find_library(METIS metis
+        "$ENV{METIS_DIR}/lib"
+        "$ENV{PARMETIS_DIR}/build/lib/libmetis"
+        "${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libmetis")
+    find_path(MFEM_INCLUDES mfem.hpp
+        "$ENV{MFEM_DIR}/include"
+        "$ENV{MFEM_DIR}"
+        "${LIBROM_DIR}/dependencies/mfem")
+    find_path(HYPRE_INCLUDES HYPRE.h
+        "$ENV{HYPRE_DIR}/include"
+        "${LIBROM_DIR}/dependencies/hypre/src/hypre/include")
+    find_path(PARMETIS_INCLUDES metis.h
+        "$ENV{PARMETIS_DIR}/metis/include"
+        "${LIBROM_DIR}/dependencies/parmetis-4.0.3/metis/include")
+    set(PYLIBROM_HAS_MFEM 1)
+endif()
 #===================== pylibROM =============================
 
 
@@ -105,30 +113,23 @@ set(CMAKE_CXX_STANDARD 14)
 
 find_package(MPI REQUIRED)
 
-set(SOURCE_DIR "bindings/pylibROM") 
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/bindings/pylibROM/pylibROM_config.h.in
+               ${CMAKE_CURRENT_SOURCE_DIR}/bindings/pylibROM/pylibROM_config.h)
+
+set(SOURCE_DIR "bindings/pylibROM")
 include_directories(
 	${SOURCE_DIR} 
 	${LIBROM_INCLUDE_DIR} 
 	${MPI_INCLUDE_PATH}
 	${MPI4PY}
 	${HDF5_C_INCLUDE_DIRS}
-	${MFEM_INCLUDES}
-	${HYPRE_INCLUDES}
-	${PARMETIS_INCLUDES}
-	${MFEM_C_INCLUDE_DIRS}
-)
-link_libraries(
-	${HDF5_LIBRARIES}
-	${MFEM}
-	${HYPRE}
-	${PARMETIS}
-	${METIS}
 )
+link_libraries(${HDF5_LIBRARIES})
 
 add_subdirectory("extern/pybind11")
 
-pybind11_add_module(_pylibROM
-    bindings/pylibROM/pylibROM.cpp
+set(PYLIBROM_SOURCES
+	bindings/pylibROM/pylibROM.cpp
 
 	bindings/pylibROM/linalg/pyMatrix.cpp
 	bindings/pylibROM/linalg/pyVector.cpp
@@ -163,14 +164,38 @@ pybind11_add_module(_pylibROM
 	bindings/pylibROM/utils/pyDatabase.hpp
 	bindings/pylibROM/utils/pyDatabase.cpp
 	bindings/pylibROM/utils/pyHDFDatabase.cpp
+	bindings/pylibROM/utils/pyHDFDatabaseMPIO.cpp
 	bindings/pylibROM/utils/pyCSVDatabase.cpp
 
-	bindings/pylibROM/mfem/pyUtilities.cpp
-	bindings/pylibROM/mfem/pyPointwiseSnapshot.cpp
-	bindings/pylibROM/mfem/pySampleMesh.cpp
-
 	bindings/pylibROM/python_utils/cpp_utils.hpp
 )
+
+if (USE_MFEM)
+    set(PYLIBROM_SOURCES ${PYLIBROM_SOURCES}
+        bindings/pylibROM/mfem/pyUtilities.cpp
+        bindings/pylibROM/mfem/pyPointwiseSnapshot.cpp
+        bindings/pylibROM/mfem/pySampleMesh.cpp)
+endif()
+
+pybind11_add_module(_pylibROM ${PYLIBROM_SOURCES})
 message("building pylibROM...")
 
+if (USE_MFEM)
+    target_include_directories(
+        _pylibROM
+        PUBLIC
+        ${MFEM_INCLUDES}
+        ${HYPRE_INCLUDES}
+        ${PARMETIS_INCLUDES}
+        ${MFEM_C_INCLUDE_DIRS})
+
+    target_link_libraries(
+        _pylibROM
+        PUBLIC
+        ${MFEM}
+        ${HYPRE}
+        ${PARMETIS}
+        ${METIS})
+endif()
+
 target_link_libraries(_pylibROM PRIVATE ROM)
diff --git a/bindings/pylibROM/__init__.py b/bindings/pylibROM/__init__.py
index 8353c0f..19641a6 100644
--- a/bindings/pylibROM/__init__.py
+++ b/bindings/pylibROM/__init__.py
@@ -6,4 +6,14 @@
 # either define/import the python routine in this file.
 # This will combine both c++ bindings/pure python routines into this module.
 
-from _pylibROM import *
+from _pylibROM.algo import *
+from _pylibROM.hyperreduction import *
+from _pylibROM.linalg import *
+
+try:
+    import _pylibROM.mfem
+    from _pylibROM.mfem import *
+except:
+    pass
+
+from _pylibROM.utils import *
diff --git a/bindings/pylibROM/linalg/pyBasisGenerator.cpp b/bindings/pylibROM/linalg/pyBasisGenerator.cpp
index be74201..30f26da 100644
--- a/bindings/pylibROM/linalg/pyBasisGenerator.cpp
+++ b/bindings/pylibROM/linalg/pyBasisGenerator.cpp
@@ -20,9 +20,9 @@ void init_BasisGenerator(pybind11::module_ &m) {
         )
         .def("isNextSample", (bool (BasisGenerator::*)(double)) &BasisGenerator::isNextSample)
         .def("updateRightSV", (bool (BasisGenerator::*)()) &BasisGenerator::updateRightSV)
-        .def("takeSample", [](BasisGenerator& self, py::array_t<double> &u_in, double time, double dt, bool add_without_increase = false) {
-            return self.takeSample(getVectorPointer(u_in), time, dt, add_without_increase);
-        }, py::arg("u_in"), py::arg("time"), py::arg("dt"), py::arg("add_without_increase") = false)
+        .def("takeSample", [](BasisGenerator& self, py::array_t<double> &u_in, bool add_without_increase = false) {
+            return self.takeSample(getVectorPointer(u_in), add_without_increase);
+        }, py::arg("u_in"), py::arg("add_without_increase") = false)
         .def("endSamples", &BasisGenerator::endSamples, py::arg("kind") = "basis")
         .def("writeSnapshot", (void (BasisGenerator::*)()) &BasisGenerator::writeSnapshot)
         .def("loadSamples", (void (BasisGenerator::*)(const std::string&, const std::string&, int, Database::formats)) &BasisGenerator::loadSamples,
@@ -39,8 +39,6 @@ void init_BasisGenerator(pybind11::module_ &m) {
         .def("getTemporalBasis", (const Matrix* (BasisGenerator::*)()) &BasisGenerator::getTemporalBasis,py::return_value_policy::reference)
         .def("getSingularValues", (const Vector* (BasisGenerator::*)()) &BasisGenerator::getSingularValues,py::return_value_policy::reference)
         .def("getSnapshotMatrix", (const Matrix* (BasisGenerator::*)()) &BasisGenerator::getSnapshotMatrix,py::return_value_policy::reference)
-        .def("getNumBasisTimeIntervals", (int (BasisGenerator::*)() const) &BasisGenerator::getNumBasisTimeIntervals)
-        .def("getBasisIntervalStartTime", (double (BasisGenerator::*)(int) const) &BasisGenerator::getBasisIntervalStartTime, py::arg("which_interval"))
         .def("getNumSamples",(int (BasisGenerator::*)() const) &BasisGenerator::getNumSamples)
         .def("__del__", [](BasisGenerator& self) { self.~BasisGenerator(); }); // Destructor
 
diff --git a/bindings/pylibROM/linalg/pyBasisReader.cpp b/bindings/pylibROM/linalg/pyBasisReader.cpp
index 4238819..0ae3042 100644
--- a/bindings/pylibROM/linalg/pyBasisReader.cpp
+++ b/bindings/pylibROM/linalg/pyBasisReader.cpp
@@ -12,54 +12,38 @@ using namespace CAROM;
 
 void init_BasisReader(pybind11::module_ &m) {
      py::class_<BasisReader>(m, "BasisReader")
-        .def(py::init<const std::string&, Database::formats>(),
+        .def(py::init<const std::string&, Database::formats, const int>(),
             py::arg("base_file_name"),
-            py::arg("file_format") = Database::formats::HDF5
+            py::arg("file_format") = Database::formats::HDF5,
+            py::arg("dim") = -1
         )
-        .def("isNewBasis",(bool (BasisReader::*)(double)) &BasisReader::isNewBasis,
-            py::arg("time"))
-        .def("getSpatialBasis",(Matrix* (BasisReader::*)(double)) &BasisReader::getSpatialBasis,
-            py::arg("time"))
-        .def("getSpatialBasis",(Matrix* (BasisReader::*)(double,int)) &BasisReader::getSpatialBasis,
-            py::arg("time"),
+        .def("getSpatialBasis",(Matrix* (BasisReader::*)()) &BasisReader::getSpatialBasis)
+        .def("getSpatialBasis",(Matrix* (BasisReader::*)(int)) &BasisReader::getSpatialBasis,
             py::arg("n"))
-        .def("getSpatialBasis",(Matrix* (BasisReader::*)(double,int,int)) &BasisReader::getSpatialBasis,
-            py::arg("time"),
+        .def("getSpatialBasis",(Matrix* (BasisReader::*)(int,int)) &BasisReader::getSpatialBasis,
             py::arg("start_col"),
             py::arg("end_col"))
-        .def("getSpatialBasis",(Matrix* (BasisReader::*)(double,double)) &BasisReader::getSpatialBasis,
-            py::arg("time"),
+        .def("getSpatialBasis",(Matrix* (BasisReader::*)(double)) &BasisReader::getSpatialBasis,
             py::arg("ef").noconvert())
-        .def("getTemporalBasis",(Matrix* (BasisReader::*)(double)) &BasisReader::getTemporalBasis,
-            py::arg("time"))
-        .def("getTemporalBasis",(Matrix* (BasisReader::*)(double,int)) &BasisReader::getTemporalBasis,
-            py::arg("time"),
+        .def("getTemporalBasis",(Matrix* (BasisReader::*)()) &BasisReader::getTemporalBasis)
+        .def("getTemporalBasis",(Matrix* (BasisReader::*)(int)) &BasisReader::getTemporalBasis,
             py::arg("n"))
-        .def("getTemporalBasis",(Matrix* (BasisReader::*)(double,int,int)) &BasisReader::getTemporalBasis,
-            py::arg("time"),
+        .def("getTemporalBasis",(Matrix* (BasisReader::*)(int,int)) &BasisReader::getTemporalBasis,
             py::arg("start_col"),
             py::arg("end_col"))
-        .def("getTemporalBasis",(Matrix* (BasisReader::*)(double,double)) &BasisReader::getTemporalBasis,
-            py::arg("time"),
+        .def("getTemporalBasis",(Matrix* (BasisReader::*)(double)) &BasisReader::getTemporalBasis,
             py::arg("ef").noconvert())
+        .def("getSingularValues",(Vector* (BasisReader::*)()) &BasisReader::getSingularValues)
         .def("getSingularValues",(Vector* (BasisReader::*)(double)) &BasisReader::getSingularValues,
-            py::arg("time"))
-        .def("getSingularValues",(Vector* (BasisReader::*)(double,double)) &BasisReader::getSingularValues,
-            py::arg("time"),
             py::arg("ef"))
-        .def("getDim", (int (BasisReader::*)(const std::string,double)) &BasisReader::getDim,
-            py::arg("kind"),
-            py::arg("time"))
-        .def("getNumSamples", (int (BasisReader::*)(const std::string,double)) &BasisReader::getNumSamples,
-            py::arg("kind"),
-            py::arg("time"))
-        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(double)) &BasisReader::getSnapshotMatrix,
-            py::arg("time"))
-        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(double,int)) &BasisReader::getSnapshotMatrix,
-            py::arg("time"),
+        .def("getDim", (int (BasisReader::*)(const std::string)) &BasisReader::getDim,
+            py::arg("kind"))
+        .def("getNumSamples", (int (BasisReader::*)(const std::string)) &BasisReader::getNumSamples,
+            py::arg("kind"))
+        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)()) &BasisReader::getSnapshotMatrix)
+        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(int)) &BasisReader::getSnapshotMatrix,
             py::arg("n"))
-        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(double,int,int)) &BasisReader::getSnapshotMatrix,
-            py::arg("time"),
+        .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(int,int)) &BasisReader::getSnapshotMatrix,
             py::arg("start_col"),
             py::arg("end_col"))
         .def("__del__", [](BasisReader& self) { self.~BasisReader(); }); // Destructor
diff --git a/bindings/pylibROM/linalg/pyMatrix.cpp b/bindings/pylibROM/linalg/pyMatrix.cpp
index 574a71f..4b03b64 100644
--- a/bindings/pylibROM/linalg/pyMatrix.cpp
+++ b/bindings/pylibROM/linalg/pyMatrix.cpp
@@ -177,7 +177,11 @@ void init_matrix(pybind11::module_ &m) {
 
         .def("transposePseudoinverse",(void (Matrix::*)()) &Matrix::transposePseudoinverse)
 
-        .def("qr_factorize",(Matrix* (Matrix::*)() const) &Matrix::qr_factorize,py::return_value_policy::take_ownership)
+        .def("qr_factorize", [](const Matrix& self) -> std::vector<std::unique_ptr<Matrix>> {
+            std::vector<std::unique_ptr<Matrix>> qr;
+            self.qr_factorize(qr);
+            return qr;
+        })
 
         // TODO (kevin): due to the difference between python and c++, technically we should not take
         //               row_pivot and row_pivot_owner as input parameters, just returning them in the end as outputs.
@@ -189,7 +193,8 @@ void init_matrix(pybind11::module_ &m) {
             return std::make_tuple(row_pivot, row_pivot_owner);
         })
 
-        .def("orthogonalize", (void (Matrix::*)()) &Matrix::orthogonalize)
+        .def("orthogonalize", (void (Matrix::*)(bool, double)) &Matrix::orthogonalize, py::arg("double_pass") = false, py::arg("zero_tol") = 1.0e-15)
+        .def("orthogonalize_last", (void (Matrix::*)(int, bool, double)) &Matrix::orthogonalize_last, py::arg("ncols") = -1, py::arg("double_pass") = false, py::arg("zero_tol") = 1.0e-15)
 
         .def("item", (const double& (Matrix::*)(int, int) const) &Matrix::item)
         .def("__getitem__", [](Matrix& self, int row, int col) { 
diff --git a/bindings/pylibROM/linalg/pyOptions.cpp b/bindings/pylibROM/linalg/pyOptions.cpp
index 302b312..2ac1424 100644
--- a/bindings/pylibROM/linalg/pyOptions.cpp
+++ b/bindings/pylibROM/linalg/pyOptions.cpp
@@ -10,10 +10,9 @@ using namespace CAROM;
 void init_Options(pybind11::module_ &m) {
 
     py::class_<Options>(m, "Options")
-        .def(py::init<int, int, int, bool, bool>(), py::arg("dim_"), py::arg("samples_per_time_interval_"),py::arg("max_time_intervals_") = -1,py::arg("update_right_SV_") = false, py::arg("write_snapshots_") = false)
+        .def(py::init<int, int, bool, bool>(), py::arg("dim_"), py::arg("max_num_samples_"),py::arg("update_right_SV_") = false, py::arg("write_snapshots_") = false)
         .def_readwrite("dim", &Options::dim)
-        .def_readwrite("samples_per_time_interval", &Options::samples_per_time_interval)
-        .def_readwrite("max_time_intervals", &Options::max_time_intervals)
+        .def_readwrite("max_num_samples", &Options::max_num_samples)
         .def_readwrite("update_right_SV", &Options::update_right_SV)
         .def_readwrite("write_snapshots", &Options::write_snapshots)
         .def_readwrite("max_basis_dimension", &Options::max_basis_dimension)
@@ -33,6 +32,7 @@ void init_Options(pybind11::module_ &m) {
         .def_readwrite("min_sampling_time_step_scale", &Options::min_sampling_time_step_scale)
         .def_readwrite("sampling_time_step_scale", &Options::sampling_time_step_scale)
         .def_readwrite("max_sampling_time_step_scale", &Options::max_sampling_time_step_scale)
+        .def_readwrite("static_svd_preserve_snapshot", &Options::static_svd_preserve_snapshot)
         .def("setMaxBasisDimension", &Options::setMaxBasisDimension, py::arg("max_basis_dimension_"))
         .def("setSingularValueTol", &Options::setSingularValueTol, py::arg("singular_value_tol_"))
         .def("setDebugMode", &Options::setDebugMode,  py::arg("debug_algorithm_"))
diff --git a/bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp b/bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp
index 275debd..8c33f21 100644
--- a/bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp
+++ b/bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp
@@ -31,8 +31,8 @@ class PyIncrementalSVD : public IncrementalSVD {
         PYBIND11_OVERRIDE(const Matrix*, IncrementalSVD, getSnapshotMatrix );
     }
 
-    bool takeSample(double* u_in, double time, bool add_without_increase) override {
-        PYBIND11_OVERLOAD(bool, IncrementalSVD, takeSample, u_in, time, add_without_increase);
+    bool takeSample(double* u_in, bool add_without_increase) override {
+        PYBIND11_OVERLOAD(bool, IncrementalSVD, takeSample, u_in, add_without_increase);
     }
 
     ~PyIncrementalSVD() override {
@@ -41,8 +41,8 @@ class PyIncrementalSVD : public IncrementalSVD {
 
 protected:
 
-    void buildInitialSVD(double* u, double time) override {
-        PYBIND11_OVERRIDE_PURE(void, IncrementalSVD, buildInitialSVD, u, time);
+    void buildInitialSVD(double* u) override {
+        PYBIND11_OVERRIDE_PURE(void, IncrementalSVD, buildInitialSVD, u);
     }
 
     void computeBasis() override {
@@ -62,18 +62,15 @@ class PyIncrementalSVD : public IncrementalSVD {
 void init_IncrementalSVD(pybind11::module_ &m) {
     py::class_<IncrementalSVD,PyIncrementalSVD>(m, "IncrementalSVD")
         .def(py::init<Options, const std::string&>())
-        .def("takeSample", [](IncrementalSVD& self, py::array_t<double> &u_in, double time,bool add_without_increase = false) {
-            bool result = self.takeSample(getVectorPointer(u_in), time, add_without_increase);
+        .def("takeSample", [](IncrementalSVD& self, py::array_t<double> &u_in, bool add_without_increase = false) {
+            bool result = self.takeSample(getVectorPointer(u_in), add_without_increase);
             return result;
-        }, py::arg("u_in"), py::arg("time"),py::arg("add_without_increase") = false)
+        }, py::arg("u_in"), py::arg("add_without_increase") = false)
         .def("getSpatialBasis", (const Matrix* (IncrementalSVD::*)()) &IncrementalSVD::getSpatialBasis)
         .def("getTemporalBasis", (const Matrix* (IncrementalSVD::*)()) &IncrementalSVD::getTemporalBasis)
         .def("getSingularValues", (const Vector* (IncrementalSVD::*)()) &IncrementalSVD::getSingularValues)
         .def("getSnapshotMatrix", (const Matrix* (IncrementalSVD::*)()) &IncrementalSVD::getSnapshotMatrix)
         .def("getDim", (int (IncrementalSVD::*)() const) &IncrementalSVD::getDim)
-        .def("getNumBasisTimeIntervals", (int (IncrementalSVD::*)() const) &IncrementalSVD::getNumBasisTimeIntervals)
-        .def("getBasisIntervalStartTime", (double (IncrementalSVD::*)(int) const) &IncrementalSVD::getBasisIntervalStartTime)
-        .def("isNewTimeInterval", (bool (IncrementalSVD::*)() const) &IncrementalSVD::isNewTimeInterval)
-        .def("increaseTimeInterval", (void (IncrementalSVD::*)()) &IncrementalSVD::increaseTimeInterval)
+        .def("getMaxNumSamples", (int (IncrementalSVD::*)() const) &IncrementalSVD::getMaxNumSamples)
         .def("getNumSamples", (int (IncrementalSVD::*)() const) &IncrementalSVD::getNumSamples);
 }
diff --git a/bindings/pylibROM/linalg/svd/pySVD.cpp b/bindings/pylibROM/linalg/svd/pySVD.cpp
index dd88a55..d27b9e8 100644
--- a/bindings/pylibROM/linalg/svd/pySVD.cpp
+++ b/bindings/pylibROM/linalg/svd/pySVD.cpp
@@ -34,8 +34,8 @@ class PySVD : public SVD {
         }
 
 
-        bool takeSample(double* u_in, double time, bool add_without_increase) override {
-        PYBIND11_OVERLOAD_PURE(bool, SVD, takeSample, u_in, time, add_without_increase);
+        bool takeSample(double* u_in, bool add_without_increase) override {
+        PYBIND11_OVERLOAD_PURE(bool, SVD, takeSample, u_in, add_without_increase);
         }
 
 
@@ -45,18 +45,15 @@ class PySVD : public SVD {
 void init_SVD(pybind11::module_ &m) {
     py::class_<SVD, PySVD>(m, "SVD")
         .def(py::init<Options>())
-        .def("takeSample", [](SVD& self, py::array_t<double> &u_in, double time,bool add_without_increase = false) {
-            return self.takeSample(getVectorPointer(u_in), time, add_without_increase);
-        }, py::arg("u_in"), py::arg("time"),py::arg("add_without_increase") = false)
+        .def("takeSample", [](SVD& self, py::array_t<double> &u_in, bool add_without_increase = false) {
+            return self.takeSample(getVectorPointer(u_in), add_without_increase);
+        }, py::arg("u_in"), py::arg("add_without_increase") = false)
         .def("getDim", (int (SVD::*)() const) &SVD::getDim)
         .def("getSpatialBasis", (const Matrix* (SVD::*)()) &SVD::getSpatialBasis)
         .def("getTemporalBasis", (const Matrix* (SVD::*)()) &SVD::getTemporalBasis)
         .def("getSingularValues", (const Vector* (SVD::*)()) &SVD::getSingularValues)
         .def("getSnapshotMatrix", (const Matrix* (SVD::*)()) &SVD::getSnapshotMatrix)
-        .def("getNumBasisTimeIntervals", (int (SVD::*)() const) &SVD::getNumBasisTimeIntervals)
-        .def("getBasisIntervalStartTime", (double (SVD::*)(int) const) &SVD::getBasisIntervalStartTime)
-        .def("isNewTimeInterval", (bool (SVD::*)() const) &SVD::isNewTimeInterval)
-        .def("increaseTimeInterval", (void (SVD::*)()) &SVD::increaseTimeInterval)
+        .def("getMaxNumSamples", (int (SVD::*)() const) &SVD::getMaxNumSamples)
         .def("getNumSamples", (int (SVD::*)() const) &SVD::getNumSamples)
         .def("__del__", [](SVD& self) {
         std::cout << "SVD instance is being destroyed" << std::endl;
diff --git a/bindings/pylibROM/linalg/svd/pyStaticSVD.cpp b/bindings/pylibROM/linalg/svd/pyStaticSVD.cpp
index daea0fe..5399fb5 100644
--- a/bindings/pylibROM/linalg/svd/pyStaticSVD.cpp
+++ b/bindings/pylibROM/linalg/svd/pyStaticSVD.cpp
@@ -14,8 +14,8 @@ class PyStaticSVD : public StaticSVD {
     using StaticSVD::StaticSVD;
 
 
-    bool takeSample(double* u_in, double time, bool add_without_increase = false) override {
-    PYBIND11_OVERRIDE(bool, StaticSVD, takeSample, u_in, time, add_without_increase);
+    bool takeSample(double* u_in, bool add_without_increase = false) override {
+    PYBIND11_OVERRIDE(bool, StaticSVD, takeSample, u_in, add_without_increase);
     }
 
 
@@ -50,19 +50,16 @@ class PyStaticSVD : public StaticSVD {
 void init_StaticSVD(pybind11::module& m) {
     py::class_<StaticSVD, PyStaticSVD>(m, "StaticSVD")
         .def(py::init(&PyStaticSVD::create), py::arg("options"))
-        .def("takeSample", [](StaticSVD& self, py::array_t<double> &u_in, double time,bool add_without_increase = false) {
-            bool result = self.takeSample(getVectorPointer(u_in), time, add_without_increase);
+        .def("takeSample", [](StaticSVD& self, py::array_t<double> &u_in, bool add_without_increase = false) {
+            bool result = self.takeSample(getVectorPointer(u_in), add_without_increase);
             return result;
-        }, py::arg("u_in"),py::arg("time"),py::arg("add_without_increase") = false)
+        }, py::arg("u_in"), py::arg("add_without_increase") = false)
         .def("getSpatialBasis", (const Matrix* (StaticSVD::*)()) &StaticSVD::getSpatialBasis,py::return_value_policy::reference_internal)
         .def("getTemporalBasis", (const Matrix* (StaticSVD::*)()) &StaticSVD::getTemporalBasis,py::return_value_policy::reference_internal)
         .def("getSingularValues", (const Vector* (StaticSVD::*)()) &StaticSVD::getSingularValues,py::return_value_policy::reference_internal)
         .def("getSnapshotMatrix", (const Matrix* (StaticSVD::*)()) &StaticSVD::getSnapshotMatrix,py::return_value_policy::reference_internal)
         .def("getDim", (int (StaticSVD::*)() const) &StaticSVD::getDim)
-        .def("getNumBasisTimeIntervals", (int (StaticSVD::*)() const) &StaticSVD::getNumBasisTimeIntervals)
-        .def("getBasisIntervalStartTime", (double (StaticSVD::*)(int) const) &StaticSVD::getBasisIntervalStartTime)
-        .def("isNewTimeInterval", (bool (StaticSVD::*)() const) &StaticSVD::isNewTimeInterval)
-        .def("increaseTimeInterval", (void (StaticSVD::*)()) &StaticSVD::increaseTimeInterval)
+        .def("getMaxNumSamples", (int (StaticSVD::*)() const) &StaticSVD::getMaxNumSamples)
         .def("getNumSamples", (int (StaticSVD::*)() const) &StaticSVD::getNumSamples);
         
 }
diff --git a/bindings/pylibROM/pylibROM.cpp b/bindings/pylibROM/pylibROM.cpp
index 6aa90ab..5f8bfdc 100644
--- a/bindings/pylibROM/pylibROM.cpp
+++ b/bindings/pylibROM/pylibROM.cpp
@@ -1,5 +1,16 @@
 #include <pybind11/pybind11.h>
 
+#include "CAROM_config.h"
+#include "pylibROM_config.h"
+
+// check that libROM has MFEM if pylibROM is using MFEM
+#ifdef PYLIBROM_HAS_MFEM
+// temporarily disabled until libROM upstream adds this option
+// #ifndef CAROM_HAS_MFEM
+// #error "libROM was not compiled with MFEM support"
+// #endif
+#endif
+
 namespace py = pybind11;
 
 //linalg
@@ -45,18 +56,22 @@ void init_Utilities(pybind11::module_ &m);
 void init_mpi_utils(pybind11::module_ &m);
 void init_Database(pybind11::module_ &m);
 void init_HDFDatabase(pybind11::module_ &m);
+void init_HDFDatabaseMPIO(pybind11::module_ &m);
 void init_CSVDatabase(pybind11::module_ &m);
 
+#ifdef PYLIBROM_HAS_MFEM
 //mfem
 void init_mfem_Utilities(pybind11::module_ &m);
 void init_mfem_PointwiseSnapshot(pybind11::module_ &m);
 void init_mfem_SampleMesh(pybind11::module_ &m);
+#endif
 
 PYBIND11_MODULE(_pylibROM, m) {
     py::module utils = m.def_submodule("utils");
     init_mpi_utils(utils);
     init_Database(utils);
     init_HDFDatabase(utils);
+    init_HDFDatabaseMPIO(utils);
     init_CSVDatabase(utils);
     
 	py::module linalg = m.def_submodule("linalg");
@@ -97,10 +112,12 @@ PYBIND11_MODULE(_pylibROM, m) {
     init_STSampling(hyperreduction);
     init_Utilities(hyperreduction);
 
+#ifdef PYLIBROM_HAS_MFEM
     py::module mfem = m.def_submodule("mfem");
     init_mfem_Utilities(mfem);
     init_mfem_PointwiseSnapshot(mfem);
     init_mfem_SampleMesh(mfem);
+#endif
 
     // py::module python_utils = m.def_submodule("python_utils");
 }
diff --git a/bindings/pylibROM/pylibROM_config.h.in b/bindings/pylibROM/pylibROM_config.h.in
new file mode 100644
index 0000000..cf1ad19
--- /dev/null
+++ b/bindings/pylibROM/pylibROM_config.h.in
@@ -0,0 +1,6 @@
+#ifndef PYLIBROM_CONFIG_H_
+#define PYLIBROM_CONFIG_H_
+
+#cmakedefine PYLIBROM_HAS_MFEM
+
+#endif
diff --git a/bindings/pylibROM/utils/mpi_utils.cpp b/bindings/pylibROM/utils/mpi_utils.cpp
index d79c52a..0c743ac 100644
--- a/bindings/pylibROM/utils/mpi_utils.cpp
+++ b/bindings/pylibROM/utils/mpi_utils.cpp
@@ -5,54 +5,14 @@
 // #include <pybind11/numpy.h>
 // #include <pybind11/operators.h>
 #include <pybind11/stl.h>
-#include "mpi.h"
-#include <mpi4py/mpi4py.h>
+
 #include "utils/mpi_utils.h"
 
+#include "mpicomm.hpp"
+
 namespace py = pybind11;
 using namespace CAROM;
 
-struct mpi4py_comm {
-    mpi4py_comm() = default;
-    mpi4py_comm(MPI_Comm value) : value(value) {}
-    operator MPI_Comm () { return value; }
-
-    MPI_Comm value;
-};
-
-namespace pybind11 {
-    namespace detail {
-        template <> struct type_caster<mpi4py_comm> {
-            public:
-            PYBIND11_TYPE_CASTER(mpi4py_comm, _("mpi4py_comm"));
-
-            // Python -> C++
-            bool load(handle src, bool) {
-                PyObject *py_src = src.ptr();
-
-                // Check that we have been passed an mpi4py communicator
-                if (PyObject_TypeCheck(py_src, &PyMPIComm_Type)) {
-                // Convert to regular MPI communicator
-                value.value = *PyMPIComm_Get(py_src);
-                } else {
-                return false;
-                }
-
-                return !PyErr_Occurred();
-            }
-
-            // C++ -> Python
-            static handle cast(mpi4py_comm src,
-                                return_value_policy /* policy */,
-                                handle /* parent */)
-            {
-                // Create an mpi4py handle
-                return PyMPIComm_New(src.value);
-            }
-        };
-    }
-} // namespace pybind11::detail
-
 void init_mpi_utils(pybind11::module_ &m) {
     // import the mpi4py API
     if (import_mpi4py() < 0) {
diff --git a/bindings/pylibROM/utils/mpicomm.hpp b/bindings/pylibROM/utils/mpicomm.hpp
new file mode 100644
index 0000000..fa4ae5a
--- /dev/null
+++ b/bindings/pylibROM/utils/mpicomm.hpp
@@ -0,0 +1,52 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "mpi.h"
+#include <mpi4py/mpi4py.h>
+
+struct mpi4py_comm {
+    mpi4py_comm() = default;
+    mpi4py_comm(MPI_Comm value) : value(value) {}
+    operator MPI_Comm () { return value; }
+
+    MPI_Comm value;
+};
+
+namespace pybind11 {
+    namespace detail {
+        template <> struct type_caster<mpi4py_comm> {
+            public:
+            PYBIND11_TYPE_CASTER(mpi4py_comm, _("mpi4py_comm"));
+
+            // Python -> C++
+            bool load(handle src, bool) {
+                PyObject *py_src = src.ptr();
+
+                // Check that we have been passed an mpi4py communicator
+                if (PyObject_TypeCheck(py_src, &PyMPIComm_Type)) {
+                // Convert to regular MPI communicator
+                value.value = *PyMPIComm_Get(py_src);
+                } else {
+                return false;
+                }
+
+                return !PyErr_Occurred();
+            }
+
+            // C++ -> Python
+            static handle cast(mpi4py_comm src,
+                                return_value_policy /* policy */,
+                                handle /* parent */)
+            {
+                if (src == MPI_COMM_NULL) {
+                    return Py_None;
+                }
+
+                // Create an mpi4py handle
+                return PyMPIComm_New(src.value);
+            }
+        };
+    }
+} // namespace pybind11::detail
\ No newline at end of file
diff --git a/bindings/pylibROM/utils/pyCSVDatabase.cpp b/bindings/pylibROM/utils/pyCSVDatabase.cpp
index 1f38b7a..b419375 100644
--- a/bindings/pylibROM/utils/pyCSVDatabase.cpp
+++ b/bindings/pylibROM/utils/pyCSVDatabase.cpp
@@ -16,34 +16,6 @@ class PyCSVDatabase : public PyDerivedDatabase<CSVDatabase> {
 public:
     using PyDerivedDatabase<CSVDatabase>::PyDerivedDatabase;
 
-    void
-    putComplexVector(
-        const std::string& file_name,
-        const std::vector<std::complex<double>>& data,
-        int nelements) override
-    {
-        PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            CSVDatabase,       /* Child class */
-            putComplexVector,        /* Name of function in C++ (must match Python name) */
-            file_name, data, nelements   /* Argument(s) */
-        );
-    }
-
-    void
-    putStringVector(
-        const std::string& file_name,
-        const std::vector<std::string>& data,
-        int nelements) override
-    {
-        PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            CSVDatabase,       /* Child class */
-            putStringVector,        /* Name of function in C++ (must match Python name) */
-            file_name, data, nelements   /* Argument(s) */
-        );
-    }
-
     // somehow this function is not virtual on c++ side. technically does not need this trampoline?
     void
     getStringVector(
@@ -81,8 +53,21 @@ void init_CSVDatabase(pybind11::module_ &m) {
     // Constructor
     csvdb.def(py::init<>());
 
-    csvdb.def("create", &CSVDatabase::create);
-    csvdb.def("open", &CSVDatabase::open);
+    csvdb.def("create", [](CSVDatabase &self, const std::string& file_name,
+                            const mpi4py_comm comm) -> bool {
+        return self.create(file_name, comm.value);
+    }, py::arg("file_name"), py::arg("comm"));
+    csvdb.def("create", [](CSVDatabase &self, const std::string& file_name) -> bool {
+        return self.create(file_name, MPI_COMM_NULL);
+    }, py::arg("file_name"));
+    csvdb.def("open", [](CSVDatabase &self, const std::string& file_name,
+                         const std::string &type, const mpi4py_comm comm) -> bool {
+        return self.open(file_name, type, comm.value);
+    }, py::arg("file_name"), py::arg("type"), py::arg("comm"));
+    csvdb.def("open", [](CSVDatabase &self, const std::string& file_name,
+                         const std::string &type) -> bool {
+        return self.open(file_name, type, MPI_COMM_NULL);
+    }, py::arg("file_name"), py::arg("type"));
     csvdb.def("close", &CSVDatabase::close);
 
     // TODO(kevin): finish binding of member functions.
diff --git a/bindings/pylibROM/utils/pyDatabase.cpp b/bindings/pylibROM/utils/pyDatabase.cpp
index b2d4d92..5ebd2dd 100644
--- a/bindings/pylibROM/utils/pyDatabase.cpp
+++ b/bindings/pylibROM/utils/pyDatabase.cpp
@@ -21,9 +21,12 @@ void init_Database(pybind11::module_ &m) {
     py::enum_<Database::formats>(db, "formats")
         .value("HDF5", Database::formats::HDF5)
         .value("CSV", Database::formats::CSV)
+        .value("HDF5_MPIO", Database::formats::HDF5_MPIO)
         .export_values();
 
-    // TODO(kevin): finish binding of member functions.
+    // https://github.com/pybind/pybind11/issues/2351
+    //db.def("create", &Database::create);
+    // // TODO(kevin): finish binding of member functions.
     db.def("getInteger", [](
         Database &self, const std::string& key)
     {
diff --git a/bindings/pylibROM/utils/pyDatabase.hpp b/bindings/pylibROM/utils/pyDatabase.hpp
index 48897e7..945b0f9 100644
--- a/bindings/pylibROM/utils/pyDatabase.hpp
+++ b/bindings/pylibROM/utils/pyDatabase.hpp
@@ -8,6 +8,7 @@
 #include <pybind11/numpy.h>
 // #include <pybind11/operators.h>
 #include <pybind11/stl.h>
+#include "utils/mpicomm.hpp"
 #include "utils/Database.h"
 
 namespace py = pybind11;
@@ -20,26 +21,30 @@ class PyDatabase : public DatabaseType
     using DatabaseType::DatabaseType; // Inherit constructors from the base class
 
     bool
-    create(const std::string& file_name) override
+    create(const std::string& file_name, const MPI_Comm comm=MPI_COMM_NULL) override
     {
         PYBIND11_OVERRIDE_PURE(
-            bool,           /* Return type */
+            bool,               /* Return type */
             DatabaseType,       /* Parent class */
-            create,         /* Name of function in C++ (must match Python name) */
-            file_name        /* Argument(s) */
+            create,             /* Name of function in C++ (must match Python name) */
+            file_name,          /* Argument(s) */
+            mpi4py_comm(comm)
         );
     }
 
     bool
     open(
         const std::string& file_name,
-        const std::string& type) override
+        const std::string& type,
+        const MPI_Comm comm=MPI_COMM_NULL) override
     {
         PYBIND11_OVERRIDE_PURE(
             bool,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            DatabaseType,          /* Parent class */
             open,                  /* Name of function in C++ (must match Python name) */
-            file_name, type         /* Argument(s) */
+            file_name,             /* Argument(s) */
+            type,
+            mpi4py_comm(comm)
         );
     }
 
@@ -48,8 +53,8 @@ class PyDatabase : public DatabaseType
     {
         PYBIND11_OVERRIDE_PURE(
             bool,                  /* Return type */
-            DatabaseType,              /* Parent class */
-            close                  /* Name of function in C++ (must match Python name) */
+            DatabaseType,          /* Parent class */
+            close,                  /* Name of function in C++ (must match Python name) */
         );
     }
 
@@ -57,13 +62,15 @@ class PyDatabase : public DatabaseType
     putIntegerArray(
         const std::string& key,
         const int* const data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
             void,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            DatabaseType,          /* Parent class */
             putIntegerArray,       /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -71,13 +78,15 @@ class PyDatabase : public DatabaseType
     putDoubleArray(
         const std::string& key,
         const double* const data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
             void,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            DatabaseType,          /* Parent class */
             putDoubleArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -85,13 +94,15 @@ class PyDatabase : public DatabaseType
     putDoubleVector(
         const std::string& key,
         const std::vector<double>& data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
-            void,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            void,                   /* Return type */
+            DatabaseType,           /* Parent class */
             putDoubleVector,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -99,13 +110,15 @@ class PyDatabase : public DatabaseType
     getIntegerArray(
         const std::string& key,
         int* data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
-            void,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            void,                   /* Return type */
+            DatabaseType,           /* Parent class */
             getIntegerArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -113,9 +126,9 @@ class PyDatabase : public DatabaseType
     getDoubleArraySize(const std::string& key) override
     {
         PYBIND11_OVERRIDE_PURE(
-            int,                  /* Return type */
-            DatabaseType,              /* Parent class */
-            getDoubleArraySize,        /* Name of function in C++ (must match Python name) */
+            int,                    /* Return type */
+            DatabaseType,           /* Parent class */
+            getDoubleArraySize,     /* Name of function in C++ (must match Python name) */
             key                     /* Argument(s) */
         );
     }
@@ -124,13 +137,15 @@ class PyDatabase : public DatabaseType
     getDoubleArray(
         const std::string& key,
         double* data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
             void,                  /* Return type */
-            DatabaseType,              /* Parent class */
+            DatabaseType,          /* Parent class */
             getDoubleArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -139,13 +154,15 @@ class PyDatabase : public DatabaseType
         const std::string& key,
         double* data,
         int nelements,
-        const std::vector<int>& idx) override
+        const std::vector<int>& idx,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
             void,                       /* Return type */
-            DatabaseType,                   /* Parent class */
+            DatabaseType,               /* Parent class */
             getDoubleArray,             /* Name of function in C++ (must match Python name) */
-            key, data, nelements, idx   /* Argument(s) */
+            key, data, nelements, idx,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -156,14 +173,16 @@ class PyDatabase : public DatabaseType
         int nelements,
         int offset,
         int block_size,
-        int stride) override
+        int stride,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE_PURE(
-            void,                       /* Return type */
-            DatabaseType,                   /* Parent class */
-            getDoubleArray,             /* Name of function in C++ (must match Python name) */
-            key, data, nelements,       /* Argument(s) */
-            offset, block_size, stride
+            void,                         /* Return type */
+            DatabaseType,                 /* Parent class */
+            getDoubleArray,               /* Name of function in C++ (must match Python name) */
+            key, data, nelements,         /* Argument(s) */
+            offset, block_size, stride,
+            distributed
         );
     }
 
@@ -175,26 +194,30 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     using PyDatabase<DerivedDatabaseType>::PyDatabase;
 
     bool
-    create(const std::string& file_name) override
+    create(const std::string& file_name,
+           const MPI_Comm comm=MPI_COMM_NULL) override
     {
         PYBIND11_OVERRIDE(
-            bool,           /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            create,         /* Name of function in C++ (must match Python name) */
-            file_name        /* Argument(s) */
+            bool,                    /* Return type */
+            DerivedDatabaseType,     /* Child class */
+            create,                  /* Name of function in C++ (must match Python name) */
+            file_name,               /* Argument(s) */
+            mpi4py_comm(comm)
         );
     }
 
     bool
     open(
         const std::string& file_name,
-        const std::string& type) override
+        const std::string& type,
+        const MPI_Comm comm=MPI_COMM_NULL) override
     {
         PYBIND11_OVERRIDE(
-            bool,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            open,                  /* Name of function in C++ (must match Python name) */
-            file_name, type         /* Argument(s) */
+            bool,                    /* Return type */
+            DerivedDatabaseType,     /* Child class */
+            open,                    /* Name of function in C++ (must match Python name) */
+            file_name, type,         /* Argument(s) */
+            mpi4py_comm(comm)
         );
     }
 
@@ -202,9 +225,9 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     close() override
     {
         PYBIND11_OVERRIDE(
-            bool,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            close                  /* Name of function in C++ (must match Python name) */
+            bool,                    /* Return type */
+            DerivedDatabaseType,     /* Child class */
+            close                    /* Name of function in C++ (must match Python name) */
         );
     }
 
@@ -212,13 +235,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     putIntegerArray(
         const std::string& key,
         const int* const data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            putIntegerArray,       /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            void,                   /* Return type */
+            DerivedDatabaseType,    /* Child class */
+            putIntegerArray,        /* Name of function in C++ (must match Python name) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -226,13 +251,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     putDoubleArray(
         const std::string& key,
         const double* const data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
             void,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
+            DerivedDatabaseType,   /* Child class */
             putDoubleArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -240,13 +267,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     putDoubleVector(
         const std::string& key,
         const std::vector<double>& data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
+            void,                   /* Return type */
+            DerivedDatabaseType,    /* Child class */
             putDoubleVector,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -254,13 +283,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     getIntegerArray(
         const std::string& key,
         int* data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
+            void,                   /* Return type */
+            DerivedDatabaseType,    /* Child class */
             getIntegerArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -268,10 +299,10 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     getDoubleArraySize(const std::string& key) override
     {
         PYBIND11_OVERRIDE(
-            int,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            getDoubleArraySize,        /* Name of function in C++ (must match Python name) */
-            key                     /* Argument(s) */
+            int,                     /* Return type */
+            DerivedDatabaseType,     /* Child class */
+            getDoubleArraySize,      /* Name of function in C++ (must match Python name) */
+            key                      /* Argument(s) */
         );
     }
 
@@ -279,13 +310,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
     getDoubleArray(
         const std::string& key,
         double* data,
-        int nelements) override
+        int nelements,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
-            void,                  /* Return type */
-            DerivedDatabaseType,       /* Child class */
-            getDoubleArray,        /* Name of function in C++ (must match Python name) */
-            key, data, nelements   /* Argument(s) */
+            void,                   /* Return type */
+            DerivedDatabaseType,    /* Child class */
+            getDoubleArray,         /* Name of function in C++ (must match Python name) */
+            key, data, nelements,   /* Argument(s) */
+            distributed
         );
     }
 
@@ -294,13 +327,15 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
         const std::string& key,
         double* data,
         int nelements,
-        const std::vector<int>& idx) override
+        const std::vector<int>& idx,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
             void,                       /* Return type */
-            DerivedDatabaseType,       /* Child class */
+            DerivedDatabaseType,        /* Child class */
             getDoubleArray,             /* Name of function in C++ (must match Python name) */
-            key, data, nelements, idx   /* Argument(s) */
+            key, data, nelements, idx,  /* Argument(s) */
+            distributed
         );
     }
 
@@ -311,14 +346,16 @@ class PyDerivedDatabase : public PyDatabase<DerivedDatabaseType> {
         int nelements,
         int offset,
         int block_size,
-        int stride) override
+        int stride,
+        const bool distributed = false) override
     {
         PYBIND11_OVERRIDE(
             void,                       /* Return type */
-            DerivedDatabaseType,       /* Child class */
+            DerivedDatabaseType,        /* Child class */
             getDoubleArray,             /* Name of function in C++ (must match Python name) */
             key, data, nelements,       /* Argument(s) */
-            offset, block_size, stride
+            offset, block_size, stride,
+            distributed
         );
     }
    
diff --git a/bindings/pylibROM/utils/pyHDFDatabase.cpp b/bindings/pylibROM/utils/pyHDFDatabase.cpp
index 769c2d4..0480eb7 100644
--- a/bindings/pylibROM/utils/pyHDFDatabase.cpp
+++ b/bindings/pylibROM/utils/pyHDFDatabase.cpp
@@ -19,8 +19,21 @@ void init_HDFDatabase(pybind11::module_ &m) {
     // Constructor
     hdfdb.def(py::init<>());
 
-    hdfdb.def("create", &HDFDatabase::create);
-    hdfdb.def("open", &HDFDatabase::open);
+    hdfdb.def("create", [](HDFDatabase &self, const std::string& file_name,
+                           const mpi4py_comm &comm) -> bool {
+        return self.create(file_name, comm.value);
+    });
+    hdfdb.def("create", [](HDFDatabase &self, const std::string& file_name) -> bool {
+        return self.create(file_name, MPI_COMM_NULL);
+    });
+    hdfdb.def("open", [](HDFDatabase &self, const std::string& file_name,
+                         const std::string &type, const mpi4py_comm &comm) -> bool {
+        return self.open(file_name, type, comm.value);
+    });
+    hdfdb.def("open", [](HDFDatabase &self, const std::string& file_name,
+                         const std::string &type) -> bool {
+        return self.open(file_name, type, MPI_COMM_NULL);
+    });
     hdfdb.def("close", &HDFDatabase::close);
 
     // TODO(kevin): finish binding of member functions.
diff --git a/bindings/pylibROM/utils/pyHDFDatabaseMPIO.cpp b/bindings/pylibROM/utils/pyHDFDatabaseMPIO.cpp
new file mode 100644
index 0000000..b10ae8c
--- /dev/null
+++ b/bindings/pylibROM/utils/pyHDFDatabaseMPIO.cpp
@@ -0,0 +1,89 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/operators.h>
+#include <pybind11/stl.h>
+#include "utils/HDFDatabaseMPIO.h"
+#include "utils/pyDatabase.hpp"
+#include "python_utils/cpp_utils.hpp"
+
+namespace py = pybind11;
+using namespace CAROM;
+
+class PyHDFDatabaseMPIO : public PyDerivedDatabase<HDFDatabaseMPIO> {
+public:
+    using PyDerivedDatabase<HDFDatabaseMPIO>::PyDerivedDatabase;
+
+};
+
+
+void init_HDFDatabaseMPIO(pybind11::module_ &m) {
+
+    py::class_<HDFDatabaseMPIO, HDFDatabase, PyDerivedDatabase<HDFDatabaseMPIO>> hdfdb(m, "HDFDatabaseMPIO");
+
+    // Constructor
+    hdfdb.def(py::init<>());
+
+    hdfdb.def("create", [](HDFDatabaseMPIO &self, const std::string& file_name,
+                           const mpi4py_comm &comm) -> bool {
+        return self.create(file_name, comm.value);
+    });
+    hdfdb.def("create", [](HDFDatabaseMPIO &self, const std::string& file_name) -> bool {
+        return self.create(file_name, MPI_COMM_NULL);
+    });
+    hdfdb.def("open", [](HDFDatabaseMPIO &self, const std::string& file_name,
+                         const std::string &type, const mpi4py_comm &comm) -> bool {
+        return self.open(file_name, type, comm.value);
+    });
+    hdfdb.def("close", &HDFDatabase::close);
+
+    // TODO(kevin): finish binding of member functions.
+    hdfdb.def("putDoubleArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, py::array_t<double> &data, int nelements)
+    {
+        self.putDoubleArray(key, getVectorPointer(data), nelements);
+    });
+    hdfdb.def("putDoubleVector", &HDFDatabase::putDoubleVector);
+
+    hdfdb.def("putInteger", &HDFDatabase::putInteger);
+    hdfdb.def("putIntegerArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, py::array_t<int> &data, int nelements)
+    {
+        self.putIntegerArray(key, getVectorPointer(data), nelements);
+    });
+
+    hdfdb.def("getIntegerArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, int nelements)
+    {
+        int *dataptr = new int[nelements];
+        self.getIntegerArray(key, dataptr, nelements);
+        return get1DArrayFromPtr(dataptr, nelements, true);
+    });
+
+    hdfdb.def("getDoubleArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, int nelements)
+    {
+        double *dataptr = new double[nelements];
+        self.getDoubleArray(key, dataptr, nelements);
+        return get1DArrayFromPtr(dataptr, nelements, true);
+    });
+
+    hdfdb.def("getDoubleArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, int nelements, const std::vector<int>& idx)
+    {
+        double *dataptr = new double[nelements];
+        self.getDoubleArray(key, dataptr, nelements, idx);
+        return get1DArrayFromPtr(dataptr, nelements, true);
+    });
+
+    hdfdb.def("getDoubleArray", [](
+        HDFDatabaseMPIO &self, const std::string& key, int nelements,
+        int offset, int block_size, int stride)
+    {
+        double *dataptr = new double[nelements];
+        self.getDoubleArray(key, dataptr, nelements, offset, block_size, stride);
+        return get1DArrayFromPtr(dataptr, nelements, true);
+    });
+
+    hdfdb.def("writeAttribute", &HDFDatabaseMPIO::writeAttribute);
+}
+
diff --git a/docker/baseline/Dockerfile b/docker/baseline/Dockerfile
index 854d14d..4b87d70 100644
--- a/docker/baseline/Dockerfile
+++ b/docker/baseline/Dockerfile
@@ -60,9 +60,9 @@ WORKDIR $LIB_DIR
 RUN git clone https://github.com/mfem/PyMFEM.git
 WORKDIR ./PyMFEM
 # --with-gslib flag does not work with later versions.
-RUN git checkout v4.5.2
+RUN git checkout v_4.7.0.1
 # We need to use this hack until we can submit a pull request to PyMFEM.
-RUN wget -O setup.py https://github.com/LLNL/pylibROM/raw/prom_parallel/extern/PyMFEM.setup.py
+#RUN wget -O setup.py https://github.com/LLNL/pylibROM/raw/prom_parallel/extern/PyMFEM.setup.py
 RUN sudo python3 setup.py install --with-parallel --with-gslib
 
 # install lldb and gdb for debugging
diff --git a/docker/librom/Dockerfile b/docker/librom/Dockerfile
index d3def6e..12f8e75 100644
--- a/docker/librom/Dockerfile
+++ b/docker/librom/Dockerfile
@@ -12,9 +12,9 @@ WORKDIR $LIB_DIR
 RUN sudo git clone https://github.com/mfem/PyMFEM.git
 WORKDIR ./PyMFEM
 # --with-gslib flag does not work with later versions.
-RUN sudo git checkout v4.5.2
+RUN sudo git checkout v_4.7.0.1
 # We need to use this hack until we can submit a pull request to PyMFEM.
-RUN sudo wget -O setup.py https://github.com/LLNL/pylibROM/raw/prom_parallel/extern/PyMFEM.setup.py
+#RUN sudo wget -O setup.py https://github.com/LLNL/pylibROM/raw/prom_parallel/extern/PyMFEM.setup.py
 RUN sudo python3 setup.py install --with-parallel --with-gslib
 
 # install libROM
@@ -22,8 +22,6 @@ WORKDIR $LIB_DIR
 RUN sudo git clone https://github.com/LLNL/libROM.git
 WORKDIR ./libROM
 RUN sudo git pull
-# pylibROM is currently based on a specific commit of the libROM.
-RUN sudo git checkout 0809d7d09dc24f0963c38fc8c0a2649948142ba0
 WORKDIR ./build
 RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=${USE_MFEM} -DMFEM_USE_GSLIB=${MFEM_USE_GSLIB}
 # RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=OFF
diff --git a/examples/prom/dg_advection_global_rom.py b/examples/prom/dg_advection_global_rom.py
index 9c016c3..ebf3e01 100644
--- a/examples/prom/dg_advection_global_rom.py
+++ b/examples/prom/dg_advection_global_rom.py
@@ -447,14 +447,14 @@ def ImplicitSolve(self, dt, x, k):
 
 # BasisGenerator for snapshot collection in offline phase 
 if offline:
-    options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots, 1,
+    options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots,
                             update_right_SV)
     generator = libROM.BasisGenerator(options, isIncremental, basisFileName)
 
 # BasisGenerator for basis construction in online phase 
 if merge:
     mergeTimer.Start()
-    options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots, 1,
+    options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots,
                             update_right_SV)
     generator = libROM.BasisGenerator(options, isIncremental, basisName)
     for paramID in range(nsets):
@@ -476,9 +476,9 @@ def ImplicitSolve(self, dt, x, k):
 if online:
     reader = libROM.BasisReader(basisName)
     if rdim != -1:
-        spatialbasis = reader.getSpatialBasis(0.0, rdim)
+        spatialbasis = reader.getSpatialBasis(rdim)
     else:
-        spatialbasis = reader.getSpatialBasis(0.0, ef)
+        spatialbasis = reader.getSpatialBasis(ef)
     numRowRB = spatialbasis.numRows()
     numColumnRB = spatialbasis.numColumns()
     if (myid == 0):
@@ -524,7 +524,7 @@ def ImplicitSolve(self, dt, x, k):
     u_centered = mfem.Vector(U.Size())
     mfem.subtract_vector(u_curr, u_init, u_centered);
     u_centered_vec = np.array((c_double * U.Size()).from_address(int(u_centered.GetData())), copy=False)
-    addSample = generator.takeSample(u_centered_vec, t, dt)
+    addSample = generator.takeSample(u_centered_vec)
 
 while not done:
     dt_real = min(dt, t_final - t)
@@ -544,7 +544,7 @@ def ImplicitSolve(self, dt, x, k):
         u_centered = mfem.Vector(U.Size())
         mfem.subtract_vector(u_curr, u_init, u_centered);
         u_centered_vec = np.array((c_double * u_centered.Size()).from_address(int(u_centered.GetData())), copy=False)
-        addSample = generator.takeSample(u_centered_vec, t, dt)
+        addSample = generator.takeSample(u_centered_vec)
 
     if done or ti % vis_steps == 0:
         if myid == 0:
diff --git a/examples/prom/dg_advection_local_rom_matrix_interp.py b/examples/prom/dg_advection_local_rom_matrix_interp.py
index 56f8b33..bbb1778 100644
--- a/examples/prom/dg_advection_local_rom_matrix_interp.py
+++ b/examples/prom/dg_advection_local_rom_matrix_interp.py
@@ -606,7 +606,7 @@ def ImplicitSolve(self, dt, x, k):
 
     # 10. Set BasisGenerator if offline
     if offline:
-        options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots, 1,
+        options = libROM.Options(fes.GetTrueVSize(), max_num_snapshots,
                             update_right_SV)
         generator = libROM.BasisGenerator(options, isIncremental, basisName)
 
@@ -614,15 +614,15 @@ def ImplicitSolve(self, dt, x, k):
         u_centered = mfem.Vector(U.Size())
         mfem.subtract_vector(u_curr, u_init, u_centered)
         u_centered_vec = np.array((c_double * U.Size()).from_address(int(u_centered.GetData())), copy=False)
-        addSample = generator.takeSample(u_centered_vec, t, dt)
+        addSample = generator.takeSample(u_centered_vec)
     
     if online:
         if not online_interp:
             reader = libROM.BasisReader(basisName)
             if rdim != -1:
-                spatialbasis = reader.getSpatialBasis(0.0, rdim)
+                spatialbasis = reader.getSpatialBasis(rdim)
             else:
-                spatialbasis = reader.getSpatialBasis(0.0, ef)
+                spatialbasis = reader.getSpatialBasis(ef)
             numRowRB = spatialbasis.numRows()
             numColumnRB = spatialbasis.numColumns()
             if (myid == 0):
@@ -684,7 +684,7 @@ def ImplicitSolve(self, dt, x, k):
                 if rdim == -1:
                     raise RuntimeError("rdim must be used for interpolation.")
 
-                parametricSpatialBasis = reader.getSpatialBasis(0.0, rdim)
+                parametricSpatialBasis = reader.getSpatialBasis(rdim)
                 numRowRB = parametricSpatialBasis.numRows()
                 numColumnRB = parametricSpatialBasis.numColumns()
                 bases.append(parametricSpatialBasis)
@@ -778,7 +778,7 @@ def ImplicitSolve(self, dt, x, k):
             u_centered = mfem.Vector(U.Size())
             mfem.subtract_vector(u_curr, u_init, u_centered)
             u_centered_vec = np.array((c_double * u_centered.Size()).from_address(int(u_centered.GetData())), copy=False)
-            addSample = generator.takeSample(u_centered_vec, t, dt)
+            addSample = generator.takeSample(u_centered_vec)
 
         if done or ti % vis_steps == 0:
             if myid == 0:
diff --git a/examples/prom/linear_elasticity_global_rom.py b/examples/prom/linear_elasticity_global_rom.py
index 7868324..8374619 100644
--- a/examples/prom/linear_elasticity_global_rom.py
+++ b/examples/prom/linear_elasticity_global_rom.py
@@ -233,13 +233,13 @@
 
 # 10. Set BasisGenerate if offline
 if offline:
-    options = la.Options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV)
+    options = la.Options(fespace.GetTrueVSize(), max_num_snapshots, update_right_SV)
     generator = la.BasisGenerator(options, isIncremental, basisFileName)
 
 # 11. The merge phase
 if merge:
     mergeTimer.Start()
-    options = la.Options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV)
+    options = la.Options(fespace.GetTrueVSize(), max_num_snapshots, update_right_SV)
     generator = la.BasisGenerator(options,isIncremental, basisName)
     for paramID in range(nsets):
         snapshot_filename = f'{basisName}{paramID}_snapshot'
@@ -338,7 +338,7 @@
     
     # 18. Take and write snapshot for ROM
     if offline:
-        addSample = generator.takeSample(X.GetDataArray(), 0.0, 0.01)
+        addSample = generator.takeSample(X.GetDataArray())
         generator.writeSnapshot()
 
 # 19. The online phase
@@ -346,7 +346,7 @@
     # 20. read the reduced basis
     assembleTimer.Start()
     reader = la.BasisReader(basisName)
-    spatialbasis = reader.getSpatialBasis(0.0)
+    spatialbasis = reader.getSpatialBasis()
     numRowRB = spatialbasis.numRows()
     numColumnRB = spatialbasis.numColumns()
     print(f'On rank {myid}, spatial basis dimension is {numRowRB} x {numColumnRB}')
diff --git a/examples/prom/nonlinear_elasticity_global_rom.py b/examples/prom/nonlinear_elasticity_global_rom.py
index f12152c..bddf692 100644
--- a/examples/prom/nonlinear_elasticity_global_rom.py
+++ b/examples/prom/nonlinear_elasticity_global_rom.py
@@ -491,7 +491,7 @@ def MergeBasis(dimFOM, nparam, max_num_snapshots, name):
     update_right_SV = False
     isIncremental = False
 
-    options = linalg.Options(dimFOM, nparam * max_num_snapshots, 1, update_right_SV)
+    options = linalg.Options(dimFOM, nparam * max_num_snapshots, update_right_SV)
     generator = linalg.BasisGenerator(options, isIncremental, "basis" + name)
 
     for paramID in range(nparam):
@@ -929,7 +929,7 @@ def run():
 
     # // 10. Create pROM object.
     if (offline):
-        options = linalg.Options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV)
+        options = linalg.Options(fespace.GetTrueVSize(), max_num_snapshots, update_right_SV)
 
         if (not x_base_only):
             basis_generator_v = linalg.BasisGenerator(options, isIncremental, basisFileName + "_V")
@@ -959,7 +959,7 @@ def run():
         else:
             readerV = linalg.BasisReader("basisV")
 
-        BV_librom = readerV.getSpatialBasis(0.0)
+        BV_librom = readerV.getSpatialBasis()
 
         if (rvdim == -1): # Change rvdim
             rvdim = BV_librom.numColumns()
@@ -972,7 +972,7 @@ def run():
             print("reduced V dim = %d\n" % rvdim)
 
         readerX = linalg.BasisReader("basisX")
-        BX_librom = readerX.getSpatialBasis(0.0)
+        BX_librom = readerX.getSpatialBasis()
 
         if (rxdim == -1): # Change rxdim
             rxdim = BX_librom.numColumns()
@@ -986,7 +986,7 @@ def run():
 
         # Hyper reduce H
         readerH = linalg.BasisReader("basisH")
-        H_librom = readerH.getSpatialBasis(0.0)
+        H_librom = readerH.getSpatialBasis()
 
         # Compute sample points
         if (hdim == -1):
@@ -1222,16 +1222,16 @@ def run():
 
             # Take samples
             if ((not x_base_only) and basis_generator_v.isNextSample(t)):
-                basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray(), t, dt)
+                basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray())
                 basis_generator_v.computeNextSampleTime(vx_diff.GetBlock(0).GetDataArray(), dvdt.GetDataArray(), t)
-                basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)
+                basis_generator_H.takeSample(oper.H_sp.GetDataArray())
 
             if (basis_generator_x.isNextSample(t)):
-                basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray(), t, dt)
+                basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray())
                 basis_generator_x.computeNextSampleTime(vx_diff.GetBlock(1).GetDataArray(), dxdt.GetDataArray(), t)
 
                 if (x_base_only):
-                    basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)
+                    basis_generator_H.takeSample(oper.H_sp.GetDataArray())
 
         if (last_step or ((ti % vis_steps) == 0)):
             if (online):
@@ -1294,15 +1294,15 @@ def run():
 
         # Take samples
         if (not x_base_only):
-            basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray(), t, dt)
+            basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray())
             basis_generator_v.writeSnapshot()
             del basis_generator_v
 
-        basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)
+        basis_generator_H.takeSample(oper.H_sp.GetDataArray())
         basis_generator_H.writeSnapshot()
         del basis_generator_H
 
-        basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray(), t, dt)
+        basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray())
         basis_generator_x.writeSnapshot()
         del basis_generator_x
 
diff --git a/examples/prom/poisson_global_rom.py b/examples/prom/poisson_global_rom.py
index 11d44ab..624ac9f 100644
--- a/examples/prom/poisson_global_rom.py
+++ b/examples/prom/poisson_global_rom.py
@@ -240,15 +240,15 @@ def run():
 
     # 10. Set BasisGenerator if offline
     if (offline):
-        options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots, 1,
-                                update_right_SV)
+        options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots,
+                                 update_right_SV)
         generator = libROM.BasisGenerator(options, isIncremental, basisFileName)
 
     # 11. The merge phase
     if (merge):
         mergeTimer.Start()
-        options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots, 1,
-                                update_right_SV)
+        options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots,
+                                 update_right_SV)
         generator = libROM.BasisGenerator(options, isIncremental, basisName)
         for paramID in range(nsets):
             snapshot_filename = "%s%d_snapshot" % (basisName, paramID)
@@ -339,7 +339,7 @@ def EvalValue(self, x):
             # NOTE: mfem Vector::GetData returns a SWIG Object of type double *.
             # To make it compatible with pybind11, we use ctypes to read data from the memory address.
             xData = np.array((c_double * X.Size()).from_address(int(X.GetData())), copy=False) # this does not copy the data.
-            addSample = generator.takeSample(xData, 0.0, 0.01)
+            addSample = generator.takeSample(xData)
             generator.writeSnapshot()
             del generator
             del options
@@ -349,7 +349,7 @@ def EvalValue(self, x):
         # 20. read the reduced basis
         assembleTimer.Start()
         reader = libROM.BasisReader(basisName)
-        spatialbasis = reader.getSpatialBasis(0.0)
+        spatialbasis = reader.getSpatialBasis()
         numRowRB = spatialbasis.numRows()
         numColumnRB = spatialbasis.numColumns()
         if (myid == 0):
diff --git a/extern/libROM b/extern/libROM
index 0809d7d..b229933 160000
--- a/extern/libROM
+++ b/extern/libROM
@@ -1 +1 @@
-Subproject commit 0809d7d09dc24f0963c38fc8c0a2649948142ba0
+Subproject commit b2299337efebf8a1affaa886c48b2172ed7e557e
diff --git a/setup.py b/setup.py
index 5d3cac1..e3f06b7 100644
--- a/setup.py
+++ b/setup.py
@@ -8,21 +8,30 @@
 # It recommends using "--config-settings", though there is no clear documentation about this.
 # For now, we enforce to use the versions < 23.3.0.
 import pkg_resources
-pkg_resources.require(['pip < 23.3.0'])
+#pkg_resources.require(['pip < 23.3.0'])
 
 from setuptools import Extension, setup, find_packages
 from setuptools.command.build_ext import build_ext
+from setuptools.command.install import install as _install
 
 # Take the global option for pre-installed librom directory.
 librom_dir = None
 install_scalapack = False
+use_mfem = True
+# TODO: fix this.. passing options through pip to setuptools with PEP517 is not working
 for arg in sys.argv:
-    if (arg[:13] == "--librom_dir="):
+    if (arg[:13] == "--librom_dir=" or arg[:13] == "--librom-dir="):
         librom_dir = arg[13:]
         sys.argv.remove(arg)
-if "--install_scalapack" in sys.argv:
-    install_scalapack = True
-    sys.argv.remove("--install_scalapack")
+    if (arg[:19] == "--install_scalapack"):
+        install_scalapack = True
+        sys.argv.remove(arg)
+    if (arg[:9] == "--no-mfem"):
+        use_mfem = False
+        sys.argv.remove(arg)
+    if (arg[:10] == "--use-mfem"):
+        use_mfem = True
+        sys.argv.remove(arg)
 
 # Convert distutils Windows platform specifiers to CMake -A arguments
 PLAT_TO_CMAKE = {
@@ -41,6 +50,38 @@ def __init__(self, name: str, sourcedir: str = "") -> None:
         self.sourcedir = os.fspath(Path(sourcedir).resolve())
 
 
+class InstallBuild(_install):
+
+    user_options = _install.user_options + [
+        ("librom-dir=", None, "Path to libROM root directory. If specified, pylibROM will use this instead of compiling its own libROM."),
+        ("install-scalapack", None, "Enables SCALAPACK when building libROM."),
+        ("use-mfem", None, "Compile pylibROM with MFEM enabled."),
+        ("no-mfem", None, "Compile pylibROM without MFEM.")
+        ]
+
+    negative_opt = dict(_install.negative_opt)
+    negative_opt.update({"no-mfem" : "use-mfem"})
+
+    def run(self):
+        _install.run(self)
+
+
+    def initialize_options(self):
+        _install.initialize_options(self)
+        self.librom_dir = None
+        self.install_scalapack = False
+        self.use_mfem = True
+
+
+    def finalize_options(self):
+        global librom_dir, install_scalapack, use_mfem
+        librom_dir = self.librom_dir
+        install_scalapack = self.install_scalapack
+        use_mfem = self.use_mfem
+
+        _install.finalize_options(self)
+
+
 class CMakeBuild(build_ext):
     def build_extension(self, ext: CMakeExtension) -> None:
         # Must be in this form due to bug in .resolve() only fixed in Python 3.10+
@@ -57,15 +98,16 @@ def build_extension(self, ext: CMakeExtension) -> None:
         # Can be set with Conda-Build, for example.
         cmake_generator = os.environ.get("CMAKE_GENERATOR", "")
 
-        global librom_dir, install_scalapack
+        global librom_dir, install_scalapack, use_mfem
         cmake_args = []
         if (librom_dir is None):
             librom_dir = os.path.dirname(os.path.realpath(__file__))
-            librom_dir += "/extern/libROM"
+            librom_dir += "/extern/libROM/"
             print("Installing libROM library: %s" % librom_dir)
             
-            librom_cmd = "cd %s && ./scripts/compile.sh -m -g -t ./cmake/toolchains/simple.cmake" % librom_dir
+            librom_cmd = "cd %s && ./scripts/compile.sh -t ./cmake/toolchains/simple.cmake" % librom_dir
             if (install_scalapack): librom_cmd += " -s"
+            if (use_mfem): librom_cmd += " -m -g"
             print("libROM installation command: %s" % librom_cmd)
             subprocess.run(
                 librom_cmd, shell=True, check=True
@@ -74,6 +116,8 @@ def build_extension(self, ext: CMakeExtension) -> None:
             print("Using pre-installed libROM library: %s" % librom_dir)
             cmake_args += [f"-DLIBROM_DIR=%s" % librom_dir]
 
+        cmake_args += ["-DUSE_MFEM={:s}".format('ON' if use_mfem else 'OFF')]
+
         # Set Python_EXECUTABLE instead if you use PYBIND11_FINDPYTHON
         # EXAMPLE_VERSION_INFO shows you how to pass a value into the C++ code
         # from Python.
@@ -152,10 +196,18 @@ def build_extension(self, ext: CMakeExtension) -> None:
             ["cmake", ext.sourcedir, *cmake_args], cwd=build_temp, check=True
         )
         subprocess.run(
-            ["cmake", "--build", ".", *build_args], cwd=build_temp, check=True
+            ["cmake", "--build", ".", "-j 8", "-v", *build_args], cwd=build_temp, check=True
         )
 
 
+    def initialize_options(self):
+        build_ext.initialize_options(self)
+
+
+    def finalize_options(self):
+        build_ext.finalize_options(self)
+
+
 # The information here can also be placed in setup.cfg - better separation of
 # logic and declaration, and simpler if you include description/version in a file.
 setup(
@@ -165,12 +217,13 @@ def build_extension(self, ext: CMakeExtension) -> None:
     # author_email="dean0x7d@gmail.com",
     description="Python Interface for LLNL libROM",
     long_description="",
-    packages=find_packages(where='bindings'),
+    packages=find_packages(where='bindings', exclude=['pylibROM.mfem'] if use_mfem == False else ['']),
     package_dir={"":"bindings"},
     # packages=['bindings/pylibROM'],
     ext_modules=[CMakeExtension("_pylibROM")],
     cmdclass={
         "build_ext": CMakeBuild,
+        "install": InstallBuild
     },
     zip_safe=False,
     extras_require={"test": ["pytest>=6.0"]},
diff --git a/tests/test_BasisGenerator.py b/tests/test_BasisGenerator.py
deleted file mode 100644
index a80d037..0000000
--- a/tests/test_BasisGenerator.py
+++ /dev/null
@@ -1,89 +0,0 @@
-import sys
-import pytest
-try:
-    # import pip-installed package
-    import pylibROM
-    import pylibROM.linalg as libROM
-    from pylibROM.utils import Database
-except ModuleNotFoundError:
-    # If pip-installed package is not found, import cmake-built package
-    sys.path.append("../build")
-    import _pylibROM as pylibROM
-    import _pylibROM.linalg as libROM
-    from _pylibROM.utils import Database
-import numpy as np 
-import h5py
-
-
-# Create an instance of BasisGenerator
-options = libROM.Options(4, 20, 3, True, True)
-incremental =   False
-basis_file_name = "basis.h5"
-file_format = Database.formats.HDF5
-generator = libROM.BasisGenerator(options, incremental, basis_file_name,file_format)
-
-
-# Test the isNextSample method
-time = 1.0  # Time of interest
-is_next = generator.isNextSample(time)
-print("isNextSample:", is_next)
-
-# Test the updateRightSV method
-update_success = generator.updateRightSV()
-print("updateRightSV success:", update_success)
-
-# Test the takeSample method
-time = 1.0  # Time of the state
-dt = 0.1  # Time step
-u_in_data = np.array([1.0, 2.0, 3.0])
-result = generator.takeSample(u_in_data, time, dt)
-print("takeSample result:", result)
-
-# Call the writeSnapshot function
-generator.writeSnapshot()
-
-# Call the loadSamples function
-# base_file_name = "basis.h5_snapshot"
-# kind = "snapshot"
-# cut_off = 10
-# db_format = Database.formats.HDF5
-# generator.loadSamples(base_file_name, kind, cut_off, db_format)
-
-# Call the computeNextSampleTime function
-u_in = [1.0, 2.0, 3.0]
-rhs_in = [0.1, 0.2, 0.3]
-time = 0.0
-next_sample_time = generator.computeNextSampleTime(u_in, rhs_in, time)
-print("computeNextSampleTime",next_sample_time)
-
-# Get the spatial basis
-spatial_basis = generator.getSpatialBasis()
-print("spatial_basis",spatial_basis.get_data())
-
-# Get the temporal basis
-temporal_basis = generator.getTemporalBasis()
-print("temporal_basis",temporal_basis.get_data())
-
-# Get the singular values
-singular_values = generator.getSingularValues()
-print("singular_values",singular_values.get_data())
-
-# Get the snapshot matrix
-snapshot_matrix = generator.getSnapshotMatrix()
-print("snapshot_matrix",snapshot_matrix.get_data())
-
-# Get the number of basis time intervals
-num_intervals = generator.getNumBasisTimeIntervals()
-print("Number of basis time intervals:", num_intervals)
-
-# Get the basis interval start time for a specific interval
-interval = 0  # Replace with the desired interval index
-interval_start_time = generator.getBasisIntervalStartTime(interval)
-print("Start time for interval", interval, ":", interval_start_time)
-
-# Get the number of samples taken
-num_samples = generator.getNumSamples()
-print("Number of samples:", num_samples)
-
-# Test the endSamples method
-generator.endSamples()
diff --git a/tests/test_Matrix.py b/tests/test_Matrix.py
deleted file mode 100644
index 7f97da4..0000000
--- a/tests/test_Matrix.py
+++ /dev/null
@@ -1,353 +0,0 @@
-import sys
-import os.path as pth
-sys.path.append(pth.join(pth.dirname(pth.abspath(__file__)), "../"))#sys.path.append("..")
-
-import build.pylibROM.linalg as libROM
-import numpy as np 
-from mpi4py import MPI 
-
-comm = MPI.COMM_WORLD
-rank = comm.Get_rank()
-size = comm.Get_size()
-print('Hello from process {} out of {}'.format(rank, size))
-
-# Create two Matrix objects
-m1 = libROM.Matrix(3,4, False,True)
-m2 = libROM.Matrix(3,4, False, False)
-
-m1.fill(2.0)
-m2.__assign__(m1)
-
-# Print the initial data for m1 and m2
-print("Initial data for m1:", m1.get_data())
-print("Initial data for m2:", m2.get_data())
-
-# Use the addition operator
-m1 += m2
-
-# Print the updated data for v1
-print("Updated data for m1 after addition:", m1.get_data())
-
-#set size 
-m1.setSize(4,5)
-
-print("Is m1 distributed:", m1.distributed())
-
-print("Is m1 balanced:", m1.balanced())
-print("Is m2 balanced:", m2.balanced())
-
-print("number of Rows in m1:", m1.numRows())
-
-print("number of columns in m1:", m1.numColumns())
-
-print("number of Distributed Rows in m1:", m1.numDistributedRows())
-
-m1.fill(3.0)
-print("Initial data for m1:", m1.get_data())
-
-# Call the getFirstNColumns method with return type Matrix*
-result1 = m1.getFirstNColumns(3)
-print("Get First 3 columns of m1",result1.get_data())  
-
-# Call the getFirstNColumns method with input parameters of n,matrix*
-m3 = libROM.Matrix()
-m1.getFirstNColumns(3, m3)
-print("Get first 3 columns of m1 to m3",m3.get_data()) 
-
-# Multiply matrix1 with matrix2
-result_matrix1 = m2.mult(m1)
-print("Matrix multiplication of m1 and m2",result_matrix1.get_data())
-
-# Multiply matrix1 with matrix2 
-result_matrix2 =  libROM.Matrix()
-m2.mult(m1,result_matrix2)
-print("The product Matrix of m1 and m2",result_matrix2.get_data())
-
-
-# Multiply matrix1 with vector2
-v1 = libROM.Vector(5, False)
-v1.fill(2.0)
-result_vector1 = m1.mult(v1)
-print("Matrix multiplication of m1 and vector v1",result_vector1.get_data())
-
-# Multiply matrix1 with matrix2 
-result_vector2 =  libROM.Vector()
-m1.mult(v1,result_vector2)
-print("The product vector of m1 and vector v1",result_vector2.get_data())
-
-# Test the first pointwise_mult function
-v2 = libROM.Vector(5, False)
-v2.fill(2.0)
-result_vector3 =  libROM.Vector(5,False)
-m1.pointwise_mult(1, v2, result_vector3)
-print("The result vector of  pointwise_multiplication of 2nd row of m1 and vector v1", result_vector3.get_data())
-
-# Test the second pointwise_mult function
-m1.pointwise_mult(1, v1)
-print("The product vector v1 after pointwise_multiplication with 2nd row of m1", v1.get_data())
-
-#Multiplies two matrices element-wise
-result_matrix=m1.elementwise_mult(m1)
-print("Result matrix of Element wise Matrix multiplication of m1 with m1",result_matrix.get_data())
-
-#Multiplies two matrices element-wise and fills result with the answer
-m3=libROM.Matrix(4,5, False,True)
-m3.fill(2.0)
-m3.elementwise_mult(m1,m3)
-print("Element wise Matrix multiplication of m3 with m1",m3.get_data())
-
-#Square every element in the matrix
-result_matrix = m1.elementwise_square()
-print("Result matrix of element wise Square multiplication of matrix m1 ",result_matrix.get_data())
-
-#Square every element in the matrix
-m1.elementwise_square(m1)
-print("Square every element in the matrix m1 ",m1.get_data())
-
-#Computes a += this*b*c
-a=libROM.Vector(4,False)
-a.fill(2.0)
-b=libROM.Vector(5,False)
-b.fill(3.0)
-c=3.0
-m3.multPlus(a,b,c)
-print("multPlus function of m3",a.get_data())
-
-# Multiplies the transpose of a Matrix with other and returns the product
-m1=m2
-result_matrix1 = m1.transposeMult(m2)
-print("Matrix multiplication of transpose of m1 and m2",result_matrix1.get_data())
-
-# Multiplies the transpose of a Matrix with other and returns the product
-result_matrix2 =  libROM.Matrix()
-m2.transposeMult(m1,result_matrix2)
-print("The product Matrix of transpose multiplication of m2 and m1",result_matrix2.get_data())
-
-
-# Multiplies the transpose of a Matrix with other vector and returns the product
-m1 = libROM.Matrix(3, 4,False,False)
-v2 = libROM.Vector(3, False)
-m1.fill(2.0)
-v2.fill(2.0)
-result_vector1 = m1.transposeMult(v2)
-print("Matrix multiplication of transpose of m1 and vector v2",result_vector1.get_data())
-
-# Multiplies the transpose of a Matrix with other vector and returns the product
-result_vector2 =  libROM.Vector()
-m1.transposeMult(v2,result_vector2) 
-print("The product vector of transpose multiplication of m1 and vector v2",result_vector2.get_data())
-
-#Computes and returns the inverse of this
-m2 = libROM.Matrix(2,2,False,False)
-m2.fill(3.0)
-m2.__setitem__(0, 0,5.0) 
-m2.__setitem__(0, 1,8.0) 
-result_matrix1 = m2.inverse()
-print("Inverse of matrix m2 (first overload):")
-print(result_matrix1.get_data())
-
-# Compute and store the inverse of m1 in the result_matrix using the second overload
-result_matrix2 = libROM.Matrix(2,2,False,False)
-m2.inverse(result_matrix2)
-print("Result matrix of inverse of matrix m2 (second overload):")
-print(result_matrix2.get_data())
-
-# Compute the inverse of m1 and store it in m1 itself using the third overload
-m2 = libROM.Matrix(2,2,False,False)
-m2.fill(3.0)
-m2.__setitem__(0, 0,5.0) 
-m2.__setitem__(0, 1,8.0) 
-m2.inverse()
-print("Matrix m2 after inverting itself (third overload):")
-print(m2.get_data())
-
-# Get a column as a Vector
-column = m1.getColumn(1)
-print("column 1 of matrix m1:", column.get_data())  
-
-#Get a column as a vector
-result_vector1=libROM.Vector()
-m1.getColumn(1,result_vector1)
-print("column 1 of matrix m1 as a vector",result_vector1.get_data())
-
-# Transpose the matrix
-print("matrix m2")
-for i in range(m2.numRows()):
-    for j in range(m2.numColumns()):
-        print(m2(i, j), end=" ")
-    print()
-m2.transpose()
-print("matrix m2 after transpose")
-# Print the transposed matrix
-for i in range(m2.numRows()):
-    for j in range(m2.numColumns()):
-        print(m2(i, j), end=" ")
-    print()
-
-# Print the transposePseudoinverse matrix
-m2.transposePseudoinverse()
-print("matrix m2 after transposePseudoinverse")
-for i in range(m2.numRows()):
-    for j in range(m2.numColumns()):
-        print(m2(i, j), end=" ")
-    print()
-
-#doubt
-# Apply qr_factorize to the matrix
-m2 = libROM.Matrix(2,2,True,False)
-m2.fill(3.0)
-m2.__setitem__(0, 0,5.0) 
-m2.__setitem__(0, 1,8.0) 
-result = m2.qr_factorize()
-
-# Print the resulting matrix
-print("qr_factorize for mtarix m2")
-for i in range(result.numRows()):
-    for j in range(result.numColumns()):
-        print(result(i, j), end=" ")
-    print()
-
-
-# Apply qrcp_pivots_transpose to the matrix
-m2 = libROM.Matrix(4,4,False,False)
-for i in range(4):
-   for j in range(4): 
-        m2.__setitem__(i,j,j) 
-print("qrcp_pivots_transpose to the matrix m2",m2.get_data())
-row_pivot = [0,0,0,0]
-row_pivot_owner = [0,0,0,0]
-row_pivots_requested = 4
-row_pivot, row_pivot_owner = m2.qrcp_pivots_transpose(row_pivot, row_pivot_owner,row_pivots_requested)
-my_rank = 0 
-for i in range(0,row_pivots_requested):
-    print("row_pivot_owner",row_pivot_owner[i])
-    print(rank)
-    assert row_pivot_owner[i] == rank
-    assert row_pivot[i] < 4
-permutation=[0,1,2,3] 
-assert np.array_equal(row_pivot, permutation)
-print("Row Pivots:", row_pivot)
-print("row_pivot_owner:", row_pivot_owner)
-
-
-# m2 = libROM.Matrix(2,2,False,False)
-# for i in range(2):
-#    for j in range(2): 
-#         m2.__setitem__(i,j,j) 
-# print("qrcp_pivots_transpose to the matrix m2",m2.get_data())
-# row_pivot = [0,0]
-# row_pivot_owner = [0,0]
-# row_pivots_requested = 2
-# row_pivot, row_pivot_owner = m2.qrcp_pivots_transpose(row_pivot, row_pivot_owner,row_pivots_requested)
-# my_rank = 0 
-# for i in range(0,row_pivots_requested):
-#     print("row_pivot_owner",i, " ",row_pivot_owner[i])
-#     print(rank)
-#     assert row_pivot_owner[i] == rank
-#     assert row_pivot[i] < 2
-# permutation=[0,1] 
-# assert np.array_equal(row_pivot, permutation)
-# print("Row Pivots:", row_pivot)
-# print("row_pivot_owner:", row_pivot_owner)
-
-
-# Apply orthogonalize to the matrix
-m2 = libROM.Matrix(2,2,False,False)
-m2.fill(3.0)
-m2.__setitem__(0, 0,5.0) 
-m2.__setitem__(0, 1,8.0) 
-m2.orthogonalize()
-print("orthogonalize to the matrix m2")
-for i in range(m2.numRows()):
-    for j in range(m2.numColumns()):
-        print(m2(i, j), end=" ")
-    print()
-
-# Set and get values using __setitem__ and __getitem__
-matrix = libROM.Matrix(3, 3,False,False)
-matrix.fill(3.0)
-matrix.__setitem__(0, 0,2.0) 
-print("Set Item (0,0) of matrix to 2.0 ",matrix.get_data())
-print("Get Item (0,0)",matrix.__getitem__(0, 0) )
-value= matrix(0, 0)
-print("value",value)
-print("call function",matrix.__call__(2,0))
-ptr=m1.getData()
-print("The storage for the Matrix's values on this processor",ptr)
-
-a=libROM.Vector(4,False)
-a.fill(2.0)
-b=libROM.Vector(5,False)
-b.fill(3.0)
-result_matrix=libROM.outerProduct(a,b)
-print("outerProduct",result_matrix.get_data())
-
-result_matrix=libROM.DiagonalMatrixFactory(a)
-print("DiagonalMatrixFactory",result_matrix.get_data())
-
-result_matrix=libROM.IdentityMatrixFactory(a)
-print("IdentityMatrixFactory",result_matrix.get_data())
-
-
-
-m1=libROM.Matrix(2,2,False,False)
-m1.__setitem__(0, 0, 4)
-m1.__setitem__(0, 1, 0)
-m1.__setitem__(1, 0, 3)
-m1.__setitem__(1, 1, -5)
-
-serialsvd1=libROM.SerialSVD(m1)
-print("U",serialsvd1.U.get_data())
-print("S",serialsvd1.S.get_data())
-print("V",serialsvd1.V.get_data())
-
-
-U=libROM.Matrix(2,2,False,False)
-S=libROM.Vector(2,False)
-V=libROM.Matrix(2,2,False,False)
-libROM.SerialSVD(m1,U,S,V)
-print("U",U.get_data())
-print("S",S.get_data())
-print("V",V.get_data())
-
-
-m1=libROM.Matrix(3,3,False,False)
-m1.__setitem__(0,0,2.0)
-m1.__setitem__(0,1,-1.0)
-m1.__setitem__(0,2,0.0)
-m1.__setitem__(1,0,-1.0)
-m1.__setitem__(1,1,3.0)
-m1.__setitem__(1,2,2.0)
-m1.__setitem__(2,0,0.0)
-m1.__setitem__(2,1,2.0)
-m1.__setitem__(2,2,4.0)
-eigenpair=libROM.SymmetricRightEigenSolve(m1)
-print("Eigen pair ev",eigenpair.ev.get_data())
-print("Eigen pair eigs",eigenpair.eigs)
-
-complexeigenpair=libROM.NonSymmetricRightEigenSolve(m1)
-print("Complex Eigen pair ev",complexeigenpair.ev_real.get_data())
-print("Complex Eigen pair ev_imaginary",complexeigenpair.ev_imaginary.get_data())
-print("Complex Eigen pair eigs",complexeigenpair.eigs)
-
-# Create the input arguments
-As = libROM.Matrix(2,2,True,False)
-As.fill(3.0)
-At = libROM.Matrix(2,2,True,False)
-At.fill(2.0)
-Bs = libROM.Matrix(2,2,True,False)
-Bs.fill(6.0)
-Bt = libROM.Matrix(2,2,True,False)
-Bt.fill(6.0)
-
-
-# Call the SpaceTimeProduct function
-m2 = libROM.Matrix(2,2,True,False)
-m2.fill(3.0)
-m1 = libROM.Matrix(2,2,True,False)
-m1.fill(5.0)
-v= libROM.Vector()
-# result = libROM.SpaceTimeProduct(As, At, Bs, Bt)
-result = libROM.SpaceTimeProduct(m1, m2, m1, m2,[1.0],False,False,False,True)
-print("SpaceTimeProduct",result.get_data())
-
diff --git a/tests/test_Options.py b/tests/test_Options.py
index a2cefd4..fa76fd2 100644
--- a/tests/test_Options.py
+++ b/tests/test_Options.py
@@ -8,16 +8,14 @@
 options = libROM.Options(4, 20)
 print("dim",options.dim)
 print("max_basis_dimension",options.max_basis_dimension)
-print("samples_per_time_interval",options.samples_per_time_interval)
-print("max_time_intervals",options.max_time_intervals)
+print("max_num_samples",options.max_num_samples)
 print("update_right_SV",options.update_right_SV)
 print("write_snapshots",options.write_snapshots)
 
 #Testing with all arguments
-options = libROM.Options(4, 20,3,True,False)
+options = libROM.Options(4, 20,True,False)
 print("dim",options.dim)
-print("samples_per_time_interval",options.samples_per_time_interval)
-print("max_time_intervals",options.max_time_intervals)
+print("max_num_samples",options.max_num_samples)
 print("update_right_SV",options.update_right_SV)
 print("write_snapshots",options.write_snapshots)
 
diff --git a/tests/test_pyBasisGenerator.py b/tests/test_pyBasisGenerator.py
index 24a29c6..a7dd75c 100644
--- a/tests/test_pyBasisGenerator.py
+++ b/tests/test_pyBasisGenerator.py
@@ -16,40 +16,38 @@
 
 
 def test_isNextSample():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     incremental = False
     basis_file_name = "basis.h5"
     file_format = Database.formats.HDF5
     generator = libROM.BasisGenerator(options, incremental, basis_file_name, file_format)
-    time = 1.0  
+    time = 1.0
     is_next = generator.isNextSample(time)
     assert is_next
 
 def test_updateRightSV():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
     update_success = generator.updateRightSV()
     assert update_success 
 
 def test_takeSample():
-    options = libROM.Options(4, 20, 3, True, True)
-    generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    time = 1.0  
-    dt = 0.1  
+    options = libROM.Options(4, 20, True, True)
+    generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5) 
     u_in_data = np.array([1.0, 2.0, 3.0])
-    result = generator.takeSample(u_in_data, time, dt)
+    result = generator.takeSample(u_in_data)
     assert result  
 
 def test_writeSnapshot():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator.takeSample(np.array([1.0, 2.0, 3.0]))
     generator.writeSnapshot()
 
 def test_computeNextSampleTime():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator1 = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator1.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator1.takeSample(np.array([1.0, 2.0, 3.0]))
     base_file_name = "test_basisgenerator_file"
     basis_writer = libROM.BasisWriter(generator1, base_file_name, Database.formats.HDF5)
     basis_writer.writeBasis("basis")
@@ -68,52 +66,38 @@ def test_computeNextSampleTime():
     generator.endSamples()
 
 def test_getSpatialBasis():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator.takeSample(np.array([1.0, 2.0, 3.0]))
     spatial_basis = generator.getSpatialBasis()
     expected_spatial_basis = np.array([[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-6.4e-323]])
     assert np.allclose(spatial_basis, expected_spatial_basis)
 
 def test_getTemporalBasis():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator.takeSample(np.array([1.0, 2.0, 3.0]))
     temporal_basis = generator.getTemporalBasis()
     expected_temporal_basis = np.array([[-1.0]])
     assert np.allclose(temporal_basis, expected_temporal_basis)
     
 def test_getSingularValues():   
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator.takeSample(np.array([1.0, 2.0, 3.0]))
     singular_values = generator.getSingularValues()
     assert(np.array_equal(singular_values.getData(),[3.7416573867739418]))
 
 def test_getSnapshotMatrix():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
+    generator.takeSample(np.array([1.0, 2.0, 3.0]))
     snapshot_matrix = generator.getSnapshotMatrix()
     expected_snapshot_basis = np.array([[1.0], [2.0], [3.0], [5.6e-322]])
     assert np.allclose(snapshot_matrix, expected_snapshot_basis)
 
-def test_getNumBasisTimeIntervals():
-    options = libROM.Options(4, 20, 3, True, True)
-    generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    num_intervals = generator.getNumBasisTimeIntervals()
-    assert num_intervals == 0
-
-def test_getBasisIntervalStartTime():
-    options = libROM.Options(4, 20, 3, True, True)
-    generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0 , 0.1)
-    interval = 0  
-    interval_start_time = generator.getBasisIntervalStartTime(interval)
-    assert interval_start_time == 1.0
-
 def test_getNumSamples():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
     num_samples = generator.getNumSamples()
     assert num_samples == 0
diff --git a/tests/test_pyBasisReader.py b/tests/test_pyBasisReader.py
index a3ee7e8..5b99b0a 100644
--- a/tests/test_pyBasisReader.py
+++ b/tests/test_pyBasisReader.py
@@ -13,9 +13,10 @@
     from _pylibROM.utils import Database
 import numpy as np 
 
-options = libROM.Options(4, 20, 3, True, True)
+options = libROM.Options(4, 20, True, True)
+options.static_svd_preserve_snapshot = True
 generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-result = generator.takeSample(np.array([1.0, 2.0, 3.0]), 0.5, 0.5)
+result = generator.takeSample(np.array([1.0, 2.0, 3.0]))
 base_file_name = "test_basisreader_file"
 basis_writer = libROM.BasisWriter(generator, base_file_name, Database.formats.HDF5)
 basis_writer.writeBasis("basis")
@@ -24,26 +25,24 @@
 generator.endSamples()
 del generator
 
-def test_isNewBasis():
-    basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
-    is_new_basis = basis_reader.isNewBasis(0.5)
-    assert is_new_basis
-    
+
 def test_getDim():
     basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
-    dim = basis_reader.getDim("basis", 0.5)
-    assert   dim == 4
-    
+    dim = basis_reader.getDim("basis")
+    assert dim == 4
+
 
 def test_getNumSamples():
     basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
-    num_samples = basis_reader.getNumSamples("basis", 0.5)
+    num_samples = basis_reader.getNumSamples("basis")
     assert num_samples == 1
-   
+
+
 def test_getSpatialBasis():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
+    options.static_svd_preserve_snapshot = True
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    result = generator.takeSample(np.array([1.0, 2.0, 3.0]), 0.5, 0.5)
+    result = generator.takeSample(np.array([1.0, 2.0, 3.0]))
     base_file_name = "test_basisreader_file"
     basis_writer = libROM.BasisWriter(generator, base_file_name, Database.formats.HDF5)
     basis_writer.writeBasis("basis")
@@ -53,60 +52,58 @@ def test_getSpatialBasis():
     del generator
     basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
 
-    spatial_basis1 = basis_reader.getSpatialBasis(0.5)
-    assert(np.allclose(spatial_basis1.getData(),[[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-4.44659081e-323]]))  
+    spatial_basis1 = basis_reader.getSpatialBasis()
+    assert(np.allclose(spatial_basis1.getData(),[[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-4.44659081e-323]]))
 
-    spatial_basis2 = basis_reader.getSpatialBasis(0.5, 1)
+    spatial_basis2 = basis_reader.getSpatialBasis(1)
     assert(np.allclose(spatial_basis2.getData(), [[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-4.44659081e-323]]))
 
-    spatial_basis3 = basis_reader.getSpatialBasis(0.5, 1, 1)
-    assert(np.allclose(spatial_basis3.getData(), [[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-4.44659081e-323]])) 
+    spatial_basis3 = basis_reader.getSpatialBasis(1, 1)
+    assert(np.allclose(spatial_basis3.getData(), [[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731], [-4.44659081e-323]]))
 
-    spatial_basis4 = basis_reader.getSpatialBasis(0.5, 0.7)
+    spatial_basis4 = basis_reader.getSpatialBasis(0.7)
     assert(np.allclose(spatial_basis4.getData(), [[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731],[-4.44659081e-323]]))
-    
+
 
 def test_getTemporalBasis():
     basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
 
-    temporal_basis1 = basis_reader.getTemporalBasis(0.5)
-    assert(np.array_equal(temporal_basis1.getData(), [[-1.0]]))  
+    temporal_basis1 = basis_reader.getTemporalBasis()
+    assert(np.array_equal(temporal_basis1.getData(), [[-1.0]]))
 
-    temporal_basis2 = basis_reader.getTemporalBasis(0.5, 1)
-    assert(np.array_equal(temporal_basis2.getData(), [[-1.0]])) 
+    temporal_basis2 = basis_reader.getTemporalBasis(1)
+    assert(np.array_equal(temporal_basis2.getData(), [[-1.0]]))
 
-    temporal_basis3 = basis_reader.getTemporalBasis(0.5, 1, 1)
+    temporal_basis3 = basis_reader.getTemporalBasis(1, 1)
     assert(np.array_equal(temporal_basis3.getData(), [[-1.0]]))
 
-    temporal_basis4 = basis_reader.getTemporalBasis(0.5, 0.7)
-    assert(np.array_equal(temporal_basis4.getData(), [[-1.0]])) 
+    temporal_basis4 = basis_reader.getTemporalBasis(0.7)
+    assert(np.array_equal(temporal_basis4.getData(), [[-1.0]]))
     
 
 def test_getSingularValues():
     basis_reader = libROM.BasisReader("test_basisreader_file", Database.formats.HDF5)
 
-    singular_values1 = basis_reader.getSingularValues(0.5)
-    assert(np.array_equal(singular_values1.getData(), [3.7416573867739418])) 
+    singular_values1 = basis_reader.getSingularValues()
+    assert(np.array_equal(singular_values1.getData(), [3.7416573867739418]))
+
+    singular_values2 = basis_reader.getSingularValues(0.7)
+    assert(np.array_equal(singular_values2.getData(), [3.7416573867739418]))
+
 
-    singular_values2 = basis_reader.getSingularValues(0.5, 0.7)
-    assert(np.array_equal(singular_values2.getData(), [3.7416573867739418]))  
-    
-    
 
 def test_getSnapshotMatrix():
     basis_reader = libROM.BasisReader("basis.h5_snapshot", Database.formats.HDF5)
 
-    snapshot_matrix1 = basis_reader.getSnapshotMatrix(0.5)
-    assert(np.allclose(snapshot_matrix1.getData(),[[-3.7416573867739418], [0.4217934441190679], [0.6326901661786019], [3.45845952e-323]])) 
-
+    snapshot_matrix1 = basis_reader.getSnapshotMatrix()
+    assert(np.allclose(snapshot_matrix1.getData(),[[1.0], [2.0], [3.0], [0.0]]))
 
-    snapshot_matrix2 = basis_reader.getSnapshotMatrix(0.5, 1)
-    assert(np.allclose(snapshot_matrix2.getData(), [[-3.7416573867739418], [0.4217934441190679], [0.6326901661786019], [3.45845952e-323]]))
+    snapshot_matrix2 = basis_reader.getSnapshotMatrix(1)
+    assert(np.allclose(snapshot_matrix2.getData(), [[1.0], [2.0], [3.0], [0.0]]))
 
+    snapshot_matrix3 = basis_reader.getSnapshotMatrix(1, 1)
+    assert(np.allclose(snapshot_matrix3.getData(), [[1.0], [2.0], [3.0], [0.0]]))
 
-    snapshot_matrix3 = basis_reader.getSnapshotMatrix(0.5, 1, 1)
-    assert(np.allclose(snapshot_matrix3.getData(), [[-3.7416573867739418], [0.4217934441190679], [0.6326901661786019], [3.45845952e-323]]))
-    
 
 if __name__ == "__main__":
     pytest.main()
diff --git a/tests/test_pyBasisWriter.py b/tests/test_pyBasisWriter.py
index 3245b51..60ab758 100644
--- a/tests/test_pyBasisWriter.py
+++ b/tests/test_pyBasisWriter.py
@@ -14,9 +14,9 @@
 import numpy as np 
 
 def test_writeBasis():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     generator = libROM.BasisGenerator(options, False, "basis.h5", Database.formats.HDF5)
-    result = generator.takeSample(np.array([1.0, 2.0, 3.0]), 1.0, 0.1) 
+    result = generator.takeSample(np.array([1.0, 2.0, 3.0])) 
     base_file_name = "basis_file"
     basis_writer = libROM.BasisWriter(generator, base_file_name, Database.formats.HDF5)
     basis_writer.writeBasis("basis")
diff --git a/tests/test_pyHDFDatabaseMPIO.py b/tests/test_pyHDFDatabaseMPIO.py
new file mode 100644
index 0000000..838d2ae
--- /dev/null
+++ b/tests/test_pyHDFDatabaseMPIO.py
@@ -0,0 +1,181 @@
+import sys
+import pytest
+import mpi4py
+try:
+    # import pip-installed package
+    import pylibROM
+    import pylibROM.linalg as libROM
+    from pylibROM.utils import Database, HDFDatabaseMPIO
+except ModuleNotFoundError:
+    # If pip-installed package is not found, import cmake-built package
+    sys.path.append("../build")
+    import _pylibROM as pylibROM
+    import _pylibROM.linalg as libROM
+    from _pylibROM.utils import Database, HDFDatabaseMPIO
+import numpy as np
+
+nrow = 123
+ncol = 21
+threshold = 1.0e-13
+
+
+def test_HDFDatabase():
+    assert mpi4py.MPI.Is_initialized()
+    comm = mpi4py.MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    nproc = comm.Get_size()
+
+    dim_rank = 2
+    nrow_local = pylibROM.split_dimension(nrow, comm)
+    global_dim, offsets = pylibROM.get_global_offsets(nrow_local, comm)
+    assert global_dim == nrow
+
+    rng = np.random.default_rng(1234)
+
+    # distribute from a global matrix to keep the same system for different nproc
+    snapshots = libROM.Matrix(nrow, ncol, False)
+    for i in range(nrow):
+        for j in range(ncol):
+            snapshots[i, j] = rng.normal(0.0, 1.0)
+    snapshots.distribute(nrow_local)
+
+    options = libROM.Options(nrow_local, ncol, True)
+    options.setMaxBasisDimension(nrow)
+    options.setRandomizedSVD(False)
+    options.setDebugMode(True)
+
+    sampler = libROM.BasisGenerator(options, False, "test_basis", Database.formats.HDF5)
+
+    sample = libROM.Vector(nrow_local, True)
+    for s in range(ncol):
+        for d in range(nrow_local):
+            sample[d] = snapshots[d, s]
+        sampler.takeSample(sample.getData())
+
+    sampler.writeSnapshot()
+    snapshot = sampler.getSnapshotMatrix()
+    sampler.endSamples()
+
+    spatial_basis = sampler.getSpatialBasis()
+
+    basis_reader = libROM.BasisReader("test_basis")
+    spatial_basis1 = basis_reader.getSpatialBasis()
+    assert spatial_basis.numRows() == spatial_basis1.numRows()
+    assert spatial_basis.numColumns() == spatial_basis1.numColumns()
+    np.testing.assert_array_almost_equal(spatial_basis.getData(), spatial_basis1.getData(), threshold)
+
+    snapshot_reader = libROM.BasisReader("test_basis_snapshot")
+    snapshot1 = snapshot_reader.getSnapshotMatrix()
+    assert snapshot.numRows() == snapshot1.numRows()
+    assert snapshot.numColumns() == snapshot1.numColumns()
+    np.testing.assert_array_almost_equal(snapshot.getData(), snapshot1.getData(), threshold)
+
+
+def test_BaseMPIOCombination():
+    assert mpi4py.MPI.Is_initialized()
+    comm = mpi4py.MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    nproc = comm.Get_size()
+
+    dim_rank = 2
+    nrow_local = pylibROM.split_dimension(nrow, comm)
+    global_dim, offsets = pylibROM.get_global_offsets(nrow_local, comm)
+    assert global_dim == nrow
+
+    base_name = "test_basis"
+    mpio_name = "test_mpio"
+
+    options = libROM.Options(nrow_local, ncol, True)
+    options.setMaxBasisDimension(nrow)
+    options.setRandomizedSVD(False)
+    options.setDebugMode(True)
+
+    sampler = libROM.BasisGenerator(options, False, mpio_name, Database.formats.HDF5_MPIO)
+
+    sampler.loadSamples(base_name + "_snapshot", "snapshot", int(1e9), Database.formats.HDF5)
+    sampler.writeSnapshot()
+    snapshot = sampler.getSnapshotMatrix()
+
+    snapshot_reader = libROM.BasisReader("test_basis_snapshot")
+    snapshot1 = snapshot_reader.getSnapshotMatrix()
+    assert snapshot.numRows() == snapshot1.numRows()
+    assert snapshot.numColumns() == snapshot1.numColumns()
+    np.testing.assert_array_almost_equal(snapshot.getData(), snapshot1.getData(), threshold)
+
+    sampler.endSamples()
+    spatial_basis = sampler.getSpatialBasis()
+
+    basis_reader = libROM.BasisReader("test_basis")
+    spatial_basis1 = basis_reader.getSpatialBasis()
+    assert spatial_basis.numRows() == spatial_basis1.numRows()
+    assert spatial_basis.numColumns() == spatial_basis1.numColumns()
+    np.testing.assert_array_almost_equal(spatial_basis.getData(), spatial_basis1.getData(), threshold)
+
+
+def test_MPIOBaseCombination():
+    assert mpi4py.MPI.Is_initialized()
+    comm = mpi4py.MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    nproc = comm.Get_size()
+
+    dim_rank = 2
+    nrow_local = pylibROM.split_dimension(nrow, comm)
+    global_dim, offsets = pylibROM.get_global_offsets(nrow_local, comm)
+    assert global_dim == nrow
+
+    base_name = "test_basis2"
+    mpio_name = "test_mpio"
+
+    options = libROM.Options(nrow_local, ncol, True)
+    options.setMaxBasisDimension(nrow)
+    options.setRandomizedSVD(False)
+    options.setDebugMode(True)
+
+    sampler = libROM.BasisGenerator(options, False, mpio_name, Database.formats.HDF5)
+
+    sampler.loadSamples(mpio_name + "_snapshot", "snapshot", int(1e9), Database.formats.HDF5_MPIO)
+    #sampler.writeSnapshot()
+    snapshot = sampler.getSnapshotMatrix()
+
+    snapshot_reader = libROM.BasisReader("test_basis_snapshot")
+    snapshot1 = snapshot_reader.getSnapshotMatrix()
+    assert snapshot.numRows() == snapshot1.numRows()
+    assert snapshot.numColumns() == snapshot1.numColumns()
+    np.testing.assert_array_almost_equal(snapshot.getData(), snapshot1.getData(), threshold)
+
+    sampler.endSamples()
+    spatial_basis = sampler.getSpatialBasis()
+
+    basis_reader = libROM.BasisReader("test_basis")
+    spatial_basis1 = basis_reader.getSpatialBasis()
+    assert spatial_basis.numRows() == spatial_basis1.numRows()
+    assert spatial_basis.numColumns() == spatial_basis1.numColumns()
+    np.testing.assert_array_almost_equal(spatial_basis.getData(), spatial_basis1.getData(), threshold)
+
+
+def test_partialGetSpatialBasis():
+    assert mpi4py.MPI.Is_initialized()
+    comm = mpi4py.MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    nproc = comm.Get_size()
+
+    nrow_local = pylibROM.split_dimension(nrow, comm)
+    global_dim, offsets = pylibROM.get_global_offsets(nrow_local, comm)
+    assert global_dim == nrow
+
+    base_name = "test_basis"
+    mpio_name = "test_mpio"
+
+    basis_reader = libROM.BasisReader(base_name)
+    spatial_basis = basis_reader.getSpatialBasis()
+
+    basis_reader1 = libROM.BasisReader(mpio_name, Database.formats.HDF5_MPIO, nrow_local)
+    spatial_basis1 = basis_reader1.getSpatialBasis()
+
+    assert spatial_basis.numRows() == spatial_basis1.numRows()
+    assert spatial_basis.numColumns() == spatial_basis1.numColumns()
+    np.testing.assert_array_almost_equal(spatial_basis.getData(), spatial_basis1.getData(), threshold)
+
+
+if __name__ == "__main__":
+    pytest.main()
diff --git a/tests/test_pyIncrementalSVD.py b/tests/test_pyIncrementalSVD.py
index 3424102..e588111 100644
--- a/tests/test_pyIncrementalSVD.py
+++ b/tests/test_pyIncrementalSVD.py
@@ -48,42 +48,12 @@ def test_getDim():
     incrementalSVD = SVD.IncrementalSVD(options, "irrelevant.txt")
     assert incrementalSVD.getDim() == 3
 
-def test_getNumBasisTimeIntervals():
+def test_getMaxNumSamples():
     options = libROM.Options(3, 4)
     options.setMaxBasisDimension(3)
     options.setIncrementalSVD(1e-1, -1.0, -1.0, -1.0)
     incrementalSVD = SVD.IncrementalSVD(options, "irrelevant.txt" )
-    num_intervals = incrementalSVD.getNumBasisTimeIntervals()
-    assert num_intervals == 0
-
-def test_getBasisIntervalStartTime():
-    options = libROM.Options(3, 4)
-    options.setMaxBasisDimension(3)
-    options.setIncrementalSVD(1e-1, -1.0, -1.0, -1.0)
-    incrementalSVD = SVD.IncrementalSVD(options, "irrelevant.txt" )
-    num_intervals = incrementalSVD.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = incrementalSVD.getBasisIntervalStartTime(interval)
-        assert interval == start_time
-
-def test_isNewTimeInterval():
-    options = libROM.Options(3, 4)
-    options.setMaxBasisDimension(3)
-    options.setIncrementalSVD(1e-1, -1.0, -1.0, -1.0)
-    incrementalSVD = SVD.IncrementalSVD(options, "irrelevant.txt" )
-    is_new_interval = incrementalSVD.isNewTimeInterval()
-    assert is_new_interval == True
-
-def test_increaseTimeInterval():
-    options = libROM.Options(3, 4)
-    options.setMaxBasisDimension(3)
-    options.setIncrementalSVD(1e-1, -1.0, -1.0, -1.0)
-    incrementalSVD = SVD.IncrementalSVD(options, "irrelevant.txt" )
-    incrementalSVD.increaseTimeInterval()
-    num_intervals = incrementalSVD.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = incrementalSVD.getBasisIntervalStartTime(interval)
-        assert interval == start_time
+    assert(incrementalSVD.getMaxNumSamples() == 4)
 
 def test_getNumSamples():
     options = libROM.Options(3, 4)
diff --git a/tests/test_pyMatrix.py b/tests/test_pyMatrix.py
index a4cf3ff..f7da0a3 100644
--- a/tests/test_pyMatrix.py
+++ b/tests/test_pyMatrix.py
@@ -467,14 +467,14 @@ def test_plus():
     result_matrix1 = m2.inverse()
     print("Inverse of matrix m2 (first overload):")
     print(result_matrix1.get_data())
-    assert result_matrix1.get_data() == [[-0.3333333333333332, 0.8888888888888886], [0.33333333333333326, -0.5555555555555554]] 
+    assert np.allclose(result_matrix1.get_data(), [[-0.3333333333333332, 0.8888888888888886], [0.33333333333333326, -0.5555555555555554]])
 
     # Compute and store the inverse of m1 in the result_matrix using the second overload
     result_matrix2 = libROM.Matrix(2,2,False,False)
     m2.inverse(result_matrix2)
     print("Result matrix of inverse of matrix m2 (second overload):")
     print(result_matrix2.get_data())
-    assert result_matrix2.get_data() == [[-0.3333333333333332, 0.8888888888888886], [0.33333333333333326, -0.5555555555555554]] 
+    assert np.allclose(result_matrix2.get_data(), [[-0.3333333333333332, 0.8888888888888886], [0.33333333333333326, -0.5555555555555554]])
 
     # Get a column as a Vector
     column = m1.getColumn(1)
@@ -527,23 +527,6 @@ def test_plus():
 
     m2 = libROM.Matrix(4,4,False,False)
 
-    # Apply orthogonalize to the matrix
-    # In parallel case, Matrix::orthogonalize() works only when matrix is distributed.
-    # This test works only for serial case, since the matrix is not distributed.
-    from mpi4py import MPI
-    if (MPI.COMM_WORLD.Get_size() == 1):
-        m2 = libROM.Matrix(2,2,False,False)
-        m2.fill(3.0)
-        m2.__setitem__(0, 0,5.0) 
-        m2.__setitem__(0, 1,8.0) 
-        m2.orthogonalize()
-        print("orthogonalize to the matrix m2")
-        for i in range(m2.numRows()):
-            for j in range(m2.numColumns()):
-                print(m2(i, j), end=" ")
-            print()
-        assert m2.get_data() == [[5.0, -0.8546160755740603],[3.0,-0.5192604003487962]]
-
     # Set and get values using __setitem__ and __getitem__
     matrix = libROM.Matrix(3, 3,False,False)
     matrix.fill(3.0)
@@ -590,9 +573,9 @@ def test_plus():
     print("U",serialsvd1.U.get_data())
     print("S",serialsvd1.S.get_data())
     print("V",serialsvd1.V.get_data())
-    assert serialsvd1.U.get_data() == [[-0.7071067811865475, 0.7071067811865475], [-0.7071067811865475, -0.7071067811865475]]
-    assert serialsvd1.S.get_data() == [6.324555320336759, 3.162277660168379]
-    assert serialsvd1.V.get_data() == [[-0.4472135954999579, -0.8944271909999159], [-0.8944271909999159, 0.4472135954999579]] 
+    np.testing.assert_allclose(serialsvd1.U.get_data(), [[-0.7071067811865475, 0.7071067811865475], [-0.7071067811865475, -0.7071067811865475]])
+    np.testing.assert_allclose(serialsvd1.S.get_data(), [6.324555320336759, 3.162277660168379])
+    np.testing.assert_allclose(serialsvd1.V.get_data(), [[-0.4472135954999579, -0.8944271909999159], [-0.8944271909999159, 0.4472135954999579]])
 
     U=libROM.Matrix(2,2,False,False)
     S=libROM.Vector(2,False)
@@ -601,9 +584,9 @@ def test_plus():
     print("U",U.get_data())
     print("S",S.get_data())
     print("V",V.get_data())
-    assert U.get_data() == [[-0.7071067811865475, 0.7071067811865475], [-0.7071067811865475, -0.7071067811865475]]
-    assert S.get_data() == [6.324555320336759, 3.162277660168379]
-    assert V.get_data() == [[-0.4472135954999579, -0.8944271909999159], [-0.8944271909999159, 0.4472135954999579]] 
+    np.testing.assert_allclose(U.get_data(), [[-0.7071067811865475, 0.7071067811865475], [-0.7071067811865475, -0.7071067811865475]])
+    np.testing.assert_allclose(S.get_data(), [6.324555320336759, 3.162277660168379])
+    np.testing.assert_allclose(V.get_data(), [[-0.4472135954999579, -0.8944271909999159], [-0.8944271909999159, 0.4472135954999579]])
 
 
     m1=libROM.Matrix(3,3,False,False)
@@ -619,16 +602,15 @@ def test_plus():
     eigenpair=libROM.SymmetricRightEigenSolve(m1)
     print("Eigen pair ev",eigenpair.ev.get_data())
     print("Eigen pair eigs",eigenpair.eigs)
-    assert eigenpair.ev.get_data() == [[0.5932333119173844, 0.7864356987513784, -0.17202653679290808], [0.6793130619863368, -0.3743619547830712, 0.6311789687764829], [-0.4319814827585531, 0.4912962635115681, 0.756320024865991]]
-    assert eigenpair.eigs == [0.8548973087995777, 2.4760236029181337, 5.669079088282289] 
-
+    np.testing.assert_allclose(eigenpair.ev.get_data(), [[0.5932333119173844, 0.7864356987513784, -0.17202653679290808], [0.6793130619863368, -0.3743619547830712, 0.6311789687764829], [-0.4319814827585531, 0.4912962635115681, 0.756320024865991]])
+    np.testing.assert_allclose(eigenpair.eigs, [0.8548973087995777, 2.4760236029181337, 5.669079088282289])
     complexeigenpair=libROM.NonSymmetricRightEigenSolve(m1)
     print("Complex Eigen pair ev",complexeigenpair.ev_real.get_data())
     print("Complex Eigen pair ev_imaginary",complexeigenpair.ev_imaginary.get_data())
     print("Complex Eigen pair eigs",complexeigenpair.eigs)
-    assert complexeigenpair.ev_real.get_data() == [[-0.5932333119173846, 0.7864356987513791, -0.17202653679290827], [-0.679313061986337, -0.37436195478307094, 0.6311789687764832], [0.43198148275855325, 0.4912962635115682, 0.7563200248659911]]
-    assert complexeigenpair.ev_imaginary.get_data() == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] 
-    assert complexeigenpair.eigs == [(0.8548973087995788+0j), (2.4760236029181346+0j), (5.669079088282289+0j)] 
+    np.testing.assert_allclose(complexeigenpair.ev_real.get_data(), [[-0.5932333119173846, 0.7864356987513791, -0.17202653679290827], [-0.679313061986337, -0.37436195478307094, 0.6311789687764832], [0.43198148275855325, 0.4912962635115682, 0.7563200248659911]])
+    np.testing.assert_allclose(complexeigenpair.ev_imaginary.get_data(), [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])
+    np.testing.assert_allclose(complexeigenpair.eigs, [(0.8548973087995788+0j), (2.4760236029181346+0j), (5.669079088282289+0j)])
 
 
     # Create the input arguments 
@@ -1078,5 +1060,98 @@ def test_distribute_and_gather():
     assert(not test.distributed())
     assert(np.array_equal(test.getData(), answer.getData()))
 
+
+def test_matrix_orthogonalize():
+    # Matrix data to orthonormalize
+    d_mat = np.array([[3.5, 7.1, 0.0, 0.0],
+                      [0.0, 1.9, 8.3, 0.0],
+                      [0.0, 0.0, 5.7, 4.6],
+                      [0.0, 0.0, 0.0, 3.2]])
+
+    # target matrix data
+    d_mat2 = np.eye(4)
+
+    matrix = libROM.Matrix(d_mat, False, False)
+    target = libROM.Matrix(d_mat2, False, False)
+
+    matrix.orthogonalize()
+
+    assert(np.allclose(matrix.getData(), target))
+
+
+def test_matrix_orthogonalize4():
+    # Matrix data to orthonormalize
+    d_mat = np.array([[3.5, 7.1, 0.0, 0.0],
+                      [0.0, 1.9, 8.3, 1.0e-14],
+                      [0.0, 0.0, 5.7, 1.0+1.0e-14],
+                      [0.0, 0.0, 0.0, 0.0]])
+
+    # target matrix data
+    d_mat2 = np.eye(4)
+    d_mat2[3][3] = 0.0
+
+    matrix = libROM.Matrix(d_mat, False, False)
+    target = libROM.Matrix(d_mat2, False, False)
+
+    matrix.orthogonalize(True)
+
+    assert(np.allclose(matrix.getData(), target))
+
+
+def test_matrix_orthogonalize_last():
+    # Matrix data to orthonormalize
+    d_mat = np.array([[1.0, 0.0, 0.0, 1.3],
+                      [0.0, 1.0, 0.0, 4.7],
+                      [0.0, 0.0, 1.0, 2.5],
+                      [0.0, 0.0, 0.0, 7.3]])
+
+    # target matrix data
+    d_mat2 = np.eye(4)
+
+    matrix = libROM.Matrix(d_mat, False, False)
+    target = libROM.Matrix(d_mat2, False, False)
+
+    matrix.orthogonalize_last()
+
+    assert(np.allclose(matrix.getData(), target))
+
+
+def test_matrix_orthogonalize_last2():
+    # Matrix data to orthonormalize
+    d_mat = np.array([[1.0, 0.0, 0.0, 1.3],
+                      [0.0, 1.0, 0.0, 4.7],
+                      [0.0, 0.0, 1.0, 2.5],
+                      [0.0, 0.0, 0.0, 7.3]])
+
+    # target matrix data
+    d_mat2 = np.eye(4)
+
+    matrix = libROM.Matrix(d_mat, False, False)
+    target = libROM.Matrix(d_mat2, False, False)
+
+    matrix.orthogonalize_last(-1, True)
+
+    assert(np.allclose(matrix.getData(), target))
+
+
+def test_matrix_orthogonalize_last4():
+    # Matrix data to orthonormalize
+    d_mat = np.array([[1.0, 0.0, 0.0, 1.3],
+                      [0.0, 1.0, 0.0, 4.7],
+                      [0.0, 0.0, 9.8, 2.5],
+                      [0.0, 0.0, 0.0, 7.3]])
+
+    # target matrix data
+    d_mat2 = np.eye(4)
+
+    matrix = libROM.Matrix(d_mat, False, False)
+    target = libROM.Matrix(d_mat2, False, False)
+
+    matrix.orthogonalize_last(3, True)
+    matrix.orthogonalize_last(4, True)
+
+    assert(np.allclose(matrix.getData(), target))
+
+
 if __name__ == '__main__':
     pytest.main()
\ No newline at end of file
diff --git a/tests/test_pyOptions.py b/tests/test_pyOptions.py
index 216c808..85c12e7 100644
--- a/tests/test_pyOptions.py
+++ b/tests/test_pyOptions.py
@@ -15,44 +15,42 @@ def test_options_constructor_two_args():
     options = libROM.Options(4, 20)
     assert options.dim == 4
     assert options.max_basis_dimension == 20
-    assert options.samples_per_time_interval == 20
-    assert options.max_time_intervals == -1
+    assert options.max_num_samples == 20
     assert not options.update_right_SV 
     assert not options.write_snapshots 
 
 def test_options_constructor_all_args():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     assert options.dim == 4
     assert options.max_basis_dimension == 20
-    assert options.samples_per_time_interval == 20
-    assert options.max_time_intervals == 3
+    assert options.max_num_samples == 20
     assert options.update_right_SV
     assert not options.write_snapshots
 
 def test_setMaxBasisDimension():
-    options = libROM.Options(4, 5, 3, True, False)
+    options = libROM.Options(4, 5, True, False)
     options.setMaxBasisDimension(10)
     assert options.max_basis_dimension == 10
 
 def test_setSingularValueTolerance():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setSingularValueTol(0.01)
     assert options.singular_value_tol == 0.01
 
 def test_setDebugMode():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setDebugMode(True)
     assert options.debug_algorithm
 
 def test_setRandomizedSVD():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setRandomizedSVD(True, randomized_subspace_dim_=5, random_seed_=42)
     assert options.randomized
     assert options.randomized_subspace_dim == 5
     assert options.random_seed == 42
 
 def test_setIncrementalSVD():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setIncrementalSVD(linearity_tol_=0.001, initial_dt_=0.1, sampling_tol_=0.001, max_time_between_samples_=1.0, fast_update_=True, skip_linearly_dependent_=False)
     assert options.linearity_tol == 0.001
     assert options.initial_dt == 0.1
@@ -62,13 +60,13 @@ def test_setIncrementalSVD():
     assert not options.skip_linearly_dependent
 
 def test_setStateIO():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setStateIO(save_state_=True, restore_state_=False)
     assert options.save_state
     assert not options.restore_state
 
 def test_setSamplingTimeStepScale():
-    options = libROM.Options(4, 20, 3, True, False)
+    options = libROM.Options(4, 20, True, False)
     options.setSamplingTimeStepScale(min_sampling_time_step_scale_=0.5, sampling_time_step_scale_=1.0, max_sampling_time_step_scale_=2.0)
     assert options.min_sampling_time_step_scale == 0.5
     assert options.sampling_time_step_scale == 1.0
diff --git a/tests/test_pyQR.py b/tests/test_pyQR.py
new file mode 100644
index 0000000..90c91db
--- /dev/null
+++ b/tests/test_pyQR.py
@@ -0,0 +1,86 @@
+import sys
+import pytest
+import mpi4py
+import numpy as np
+try:
+    # import pip-installed package
+    import pylibROM
+    import pylibROM.linalg as libROM
+except ModuleNotFoundError:
+    # If pip-installed package is not found, import cmake-built package
+    sys.path.append("../build")
+    import _pylibROM as pylibROM
+    import _pylibROM.linalg as libROM
+
+
+def test_MatrixQR():
+    assert mpi4py.MPI.Is_initialized()
+    comm = mpi4py.MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    num_procs = comm.Get_size()
+
+    num_total_rows = 5
+    num_columns = 3
+    loc_num_rows = pylibROM.split_dimension(num_total_rows, comm)
+    total_rows, row_offset = pylibROM.get_global_offsets(loc_num_rows, comm)
+    assert total_rows == num_total_rows
+
+    q_data = [
+        3.08158946098238906153E-01,    -9.49897947980619661301E-02,    -4.50691774108525788911E-01,
+        -1.43697905723455976457E-01,    9.53289043424090820622E-01,    8.77767692937209131898E-02,
+        -2.23655845793717528158E-02,    -2.10628953513210204207E-01,   8.42235962392685943989E-01,
+        -7.29903965154318323805E-01,    -1.90917141788945754488E-01,   -2.77280930877637610266E-01,
+        -5.92561353877168350834E-01,    -3.74570084880578441089E-02,   5.40928141934190823137E-02
+    ]
+
+    r_data = [
+        -1.78651649346571794741E-01,     5.44387957786310106023E-01,    -8.19588518467042281834E-01,
+        0.0,                             -3.13100149275943651084E-01,   -9.50441422536040881122E-04,
+        0.0,                             0.0,                            5.72951792961765460355E-01
+    ]
+
+    exactQ = libROM.Matrix(loc_num_rows, num_columns, True)
+    exactR = libROM.Matrix(np.asarray(r_data).reshape((3,3)), False, True)
+
+    for i in range(loc_num_rows):
+        for j in range(num_columns):
+            exactQ[i,j] = q_data[((i + row_offset[rank]) * num_columns) + j]
+
+    assert exactQ.numRows() == loc_num_rows
+    assert exactQ.numColumns() == num_columns
+
+    # Verify that the columns of Q are orthonormal
+    id = exactQ.transposeMult(exactQ)
+    assert id.numRows() == num_columns
+    assert id.numColumns() == num_columns
+
+    maxError = np.max(np.abs(id.getData() - np.eye(num_columns)))
+    np.testing.assert_almost_equal(maxError, 0.0, 15)
+
+
+    # Compute A = QR
+    # A = exactQ.mult(exactR) # need PR28 fix for this syntax
+    A = libROM.Matrix(num_columns, num_columns, True)
+    exactQ.mult(exactR, A)
+
+    # Factorize A
+    QRfactors = A.qr_factorize()
+    assert len(QRfactors) == 2
+
+    Q = QRfactors[0]
+    R = QRfactors[1]
+    assert Q.numRows() == loc_num_rows
+    assert Q.numColumns() == num_columns
+    assert R.numRows() == num_columns
+    assert R.numColumns() == num_columns
+
+    # Verify that Q == -exactQ and R == -exactR
+    maxError = np.max(np.abs(exactQ.getData() + Q.getData()))
+    np.testing.assert_almost_equal(maxError, 0.0, 15)
+
+    maxError = np.max(np.abs(exactR.getData() + R.getData()))
+    np.testing.assert_almost_equal(maxError, 0.0, 15)
+
+
+if __name__ == '__main__':
+    pytest.main()
diff --git a/tests/test_pySTSampling.py b/tests/test_pySTSampling.py
index 88139be..cb9d3b6 100644
--- a/tests/test_pySTSampling.py
+++ b/tests/test_pySTSampling.py
@@ -5,6 +5,8 @@
 import numpy as np 
 import _pylibROM.hyperreduction as hyperreduction
 
+
+@pytest.mark.skip("TODO")
 def test_SpaceTimeSampling():
         # Prepare test data
         num_f_basis_vectors_used = 5
@@ -12,19 +14,23 @@ def test_SpaceTimeSampling():
         num_rows = 10
         num_t_samples_req=5
         num_s_samples_req= 10
-        s_basis = linalg.Matrix(num_rows,num_cols,False,False)
-        t_basis = linalg.Matrix(num_rows,num_cols,True,False)
+        s_basis = linalg.Matrix(num_rows,num_cols,True,False)
+        t_basis = linalg.Matrix(num_rows,num_cols,False,False)
         f_sampled_row = np.zeros(num_s_samples_req, dtype=np.int32)
         f_sampled_rows_per_proc = np.zeros(1, dtype=np.int32)
         s_basis_sampled = linalg.Matrix(num_s_samples_req,num_t_samples_req,False,False)
         myid = 0
         num_procs = 1
 
-        t_samples= hyperreduction.SpaceTimeSampling(s_basis, t_basis, num_f_basis_vectors_used, f_sampled_row,f_sampled_rows_per_proc, s_basis_sampled, myid, num_procs,num_t_samples_req,num_s_samples_req,False)
+        t_samples= hyperreduction.SpaceTimeSampling(s_basis, t_basis, num_f_basis_vectors_used,
+                                                    f_sampled_row,f_sampled_rows_per_proc, s_basis_sampled,
+                                                    myid, num_procs,num_t_samples_req,num_s_samples_req,False)
         assert np.all(t_samples == [0, 1, 2, 3, 4])
         assert np.all(f_sampled_row == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
         assert np.all(f_sampled_rows_per_proc == [10])
 
+
+@pytest.mark.skip("TODO")
 def test_GetSampledSpaceTimeBasis():
         t_samples = [1,2,3]
         t_basis = linalg.Matrix(3,5,True,False)
@@ -34,7 +40,8 @@ def test_GetSampledSpaceTimeBasis():
         f_basis_sampled_inv = linalg.Matrix(6, 2, False,False)
         hyperreduction.GetSampledSpaceTimeBasis(t_samples, t_basis, s_basis_sampled, f_basis_sampled_inv)
         assert np.allclose(f_basis_sampled_inv.getData(),[[1728.0, 6.0], [1728.0, 6.0], [4.751093e-318, 3.297888e-320], [1728.0, 6.0], [1728.0, 6.0], [4.751093e-318, 3.297888e-320]])
-     
+
+
 if __name__ == '__main__':
     pytest.main()
 
diff --git a/tests/test_pySVD.py b/tests/test_pySVD.py
index 7eb7be4..d98de44 100644
--- a/tests/test_pySVD.py
+++ b/tests/test_pySVD.py
@@ -14,42 +14,18 @@
 import numpy as np 
 
 def test_getDim():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     svd = SVD.SVD(options)
     dim = svd.getDim()
     assert dim == 4
 
-def test_getNumBasisTimeIntervals():
-    options = libROM.Options(4, 20, 3, True, True)
+def test_getMaxNumSamples():
+    options = libROM.Options(4, 20, True, True)
     svd = SVD.SVD(options)
-    num_intervals = svd.getNumBasisTimeIntervals()
-    assert num_intervals == 0
-
-def test_getBasisIntervalStartTime():
-    options = libROM.Options(4, 20, 3, True, True)
-    svd = SVD.SVD(options)
-    num_intervals = svd.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = svd.getBasisIntervalStartTime(interval)
-        assert interval == start_time
-
-def test_isNewTimeInterval():
-    options = libROM.Options(4, 20, 3, True, True)
-    svd = SVD.SVD(options)
-    is_new_interval = svd.isNewTimeInterval()
-    assert is_new_interval == True
-
-def test_increaseTimeInterval():
-    options = libROM.Options(4, 20, 3, True, True)
-    svd = SVD.SVD(options)
-    svd.increaseTimeInterval()
-    num_intervals = svd.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = svd.getBasisIntervalStartTime(interval)
-        assert interval == start_time
+    assert(svd.getMaxNumSamples() == 20)
 
 def test_getNumSamples():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     svd = SVD.SVD(options)
     num_samples = svd.getNumSamples()
     assert num_samples == 0
diff --git a/tests/test_pyStaticSVD.py b/tests/test_pyStaticSVD.py
index 24176ad..4428b2f 100644
--- a/tests/test_pyStaticSVD.py
+++ b/tests/test_pyStaticSVD.py
@@ -14,96 +14,67 @@
 import numpy as np 
 
 def test_getDim():
-    options = libROM.Options(3, 10, 3, True, True)
+    options = libROM.Options(3, 10, True, True)
     staticsvd=SVD.StaticSVD(options)
     dim = staticsvd.getDim()
     assert dim == 3
 
-def test_getNumBasisTimeIntervals():
-    options = libROM.Options(3, 10, 3, True, True)
-    staticsvd=SVD.StaticSVD(options)
-    num_intervals = staticsvd.getNumBasisTimeIntervals()
-    assert num_intervals == 0
-
-def test_getBasisIntervalStartTime():
-    options = libROM.Options(3, 10, 3, True, True)
-    staticsvd=SVD.StaticSVD(options)
-    num_intervals = staticsvd.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = staticsvd.getBasisIntervalStartTime(interval)
-        assert interval == start_time
-
-def test_isNewTimeInterval():
-    options = libROM.Options(3, 10, 3, True, True)
-    staticsvd=SVD.StaticSVD(options)
-    is_new_interval = staticsvd.isNewTimeInterval()
-    assert is_new_interval == True
-
-def test_increaseTimeInterval():
-    options = libROM.Options(3, 10, 3, True, True)
-    staticsvd=SVD.StaticSVD(options)
-    staticsvd.increaseTimeInterval()
-    num_intervals = staticsvd.getNumBasisTimeIntervals()
-    for interval in range(num_intervals):
-        start_time = staticsvd.getBasisIntervalStartTime(interval)
-        assert interval == start_time
-
 def test_getNumSamples():
-    options = libROM.Options(4, 20, 3, True, True)
+    options = libROM.Options(4, 20, True, True)
     staticsvd = SVD.StaticSVD(options)
     num_samples = staticsvd.getNumSamples()
     assert num_samples == 0
 
+def test_getMaxNumSamples():
+    options = libROM.Options(4, 20, True, True)
+    staticsvd = SVD.StaticSVD(options)
+    assert(staticsvd.getMaxNumSamples() == 20)
+
 def test_takeSample():
-    options = libROM.Options(3, 10, 3, True, True)
+    options = libROM.Options(3, 10, True, True)
     staticsvd = SVD.StaticSVD(options)
     u_in = np.array([1.0, 2.0, 3.0, 4.0])
-    time = 0.5
     add_without_increase = False
-    result = staticsvd.takeSample(u_in, time, add_without_increase)
+    result = staticsvd.takeSample(u_in, add_without_increase)
     assert result == True
 
 def test_getSpatialBasis():
-    options = libROM.Options(3, 10, 3, True, True)
+    options = libROM.Options(3, 10, True, True)
     staticsvd = SVD.StaticSVD(options)
     u_in = np.array([1.0, 2.0, 3.0, 4.0])
-    time = 0.5
     add_without_increase = False
-    staticsvd.takeSample(u_in, time, add_without_increase)
+    staticsvd.takeSample(u_in, add_without_increase)
     spatial_basis = staticsvd.getSpatialBasis()
     assert(np.array_equal(spatial_basis.getData(), [[-0.2672612419124243], [-0.5345224838248487], [-0.8017837257372731]]))
 
 def test_getTemporalBasis():
-    options = libROM.Options(3, 10, 3, True, True)
+    options = libROM.Options(3, 10, True, True)
     staticsvd = SVD.StaticSVD(options)
     u_in = np.array([1.0, 2.0, 3.0, 4.0])
-    time = 0.5
     add_without_increase = False
-    staticsvd.takeSample(u_in, time, add_without_increase)
+    staticsvd.takeSample(u_in, add_without_increase)
     temporal_basis = staticsvd.getTemporalBasis()
     assert(np.array_equal(temporal_basis.getData(), [[-1.0]]))
 
 def test_getSingularValues():
-    options = libROM.Options(3, 10, 3, True, True)
+    options = libROM.Options(3, 10, True, True)
     staticsvd = SVD.StaticSVD(options)
     u_in = np.array([1.0, 2.0, 3.0, 4.0])
-    time = 0.5
     add_without_increase = False
-    staticsvd.takeSample(u_in, time, add_without_increase)
+    staticsvd.takeSample(u_in, add_without_increase)
     singular_values = staticsvd.getSingularValues()
     assert(np.array_equal(singular_values.getData(), [3.7416573867739418]))
 
 def test_getSnapshotMatrix():
-    options = libROM.Options(5, 10, 5, True, True)
+    options = libROM.Options(5, 10, True, True)
     staticsvd = SVD.StaticSVD(options)
     input_snapshot = np.array([[0.5377, 1.8339, -2.2588, 0.8622, 0.3188],
                            [-1.3077, -0.4336, 0.3426, 3.5784, 2.7694],
                            [-1.3499, 3.0349, 0.7254, -0.0631, 0.7147]])
     input_snapshot = input_snapshot.T
-    time = 0.5
     add_without_increase = False
     for k in range(input_snapshot.shape[1]):
-         result = staticsvd.takeSample(input_snapshot[:,k], time, add_without_increase)
+         result = staticsvd.takeSample(input_snapshot[:,k], add_without_increase)
     snapshot_matrix = staticsvd.getSnapshotMatrix()
     assert(np.array_equal(snapshot_matrix, input_snapshot))