diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp
index d0b6514d3..4918b1ba1 100644
--- a/lib/linalg/Matrix.cpp
+++ b/lib/linalg/Matrix.cpp
@@ -1788,47 +1788,88 @@ const
 }
 
 void
-Matrix::orthogonalize()
+Matrix::orthogonalize(double zero_tol)
 {
-    for (int work = 1; work < d_num_cols; ++work) {
-        double tmp;
-        for (int col = 0; col < work; ++col) {
+    for (int work = 0; work < d_num_cols; ++work)
+    {
+        // Orthogonalize the column.
+        for (int col = 0; col < work; ++col)
+        {
             double factor = 0.0;
-            tmp = 0.0;
-            for (int i = 0; i < d_num_rows; ++i) {
-                tmp += item(i, col)*item(i, work);
-            }
-            if (d_num_procs > 1) {
-                MPI_Allreduce(&tmp,
-                              &factor,
-                              1,
-                              MPI_DOUBLE,
-                              MPI_SUM,
-                              MPI_COMM_WORLD);
-            }
-            else {
-                factor = tmp;
-            }
 
-            for (int i = 0; i < d_num_rows; ++i) {
-                item(i, work) -= factor*item(i, col);
+            for (int i = 0; i < d_num_rows; ++i)
+                factor += item(i, col) * item(i, work);
+
+            if (d_distributed && d_num_procs > 1)
+            {
+                CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1,
+                                            MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
             }
+            for (int i = 0; i < d_num_rows; ++i)
+                item(i, work) -= factor * item(i, col);
         }
+
+        // Normalize the column.
         double norm = 0.0;
-        tmp = 0.0;
-        for (int i = 0; i < d_num_rows; ++i) {
-            tmp += item(i, work)*item(i, work);
-        }
-        if (d_num_procs > 1) {
-            MPI_Allreduce(&tmp, &norm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+        for (int i = 0; i < d_num_rows; ++i)
+            norm += item(i, work) * item(i, work);
+
+        if (d_distributed && d_num_procs > 1)
+        {
+            CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE,
+                                        MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
         }
-        else {
-            norm = tmp;
+        if (norm > zero_tol)
+        {
+            norm = 1.0 / sqrt(norm);
+            for (int i = 0; i < d_num_rows; ++i)
+                item(i, work) *= norm;
         }
-        norm = sqrt(norm);
-        for (int i = 0; i < d_num_rows; ++i) {
-            item(i, work) /= norm;
+    }
+}
+
+void
+Matrix::orthogonalize_last(int ncols, double zero_tol)
+{
+    if (ncols == -1) ncols = d_num_cols;
+    CAROM_VERIFY((ncols > 0) && (ncols <= d_num_cols));
+
+    const int last_col = ncols - 1; // index of column to be orthonormalized
+
+    // Orthogonalize the column.
+    for (int col = 0; col < last_col; ++col)
+    {
+        double factor = 0.0;
+
+        for (int i = 0; i < d_num_rows; ++i)
+            factor += item(i, col) * item(i, last_col);
+
+        if (d_distributed && d_num_procs > 1)
+        {
+            CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &factor, 1, MPI_DOUBLE,
+                                        MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
         }
+        for (int i = 0; i < d_num_rows; ++i)
+            item(i, last_col) -= factor * item(i, col);
+    }
+
+    // Normalize the column.
+    double norm = 0.0;
+
+    for (int i = 0; i < d_num_rows; ++i)
+        norm += item(i, last_col) * item(i, last_col);
+
+    if (d_distributed && d_num_procs > 1)
+    {
+        CAROM_VERIFY( MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE,
+                                    MPI_SUM, MPI_COMM_WORLD) == MPI_SUCCESS );
+    }
+    if (norm > zero_tol)
+    {
+        norm = 1.0 / sqrt(norm);
+        for (int i = 0; i < d_num_rows; ++i)
+            item(i, last_col) *= norm;
     }
 }
 
diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h
index 4f7f74aba..e5f5bc988 100644
--- a/lib/linalg/Matrix.h
+++ b/lib/linalg/Matrix.h
@@ -941,10 +941,33 @@ class Matrix
                           int  pivots_requested) const;
 
     /**
-     * @brief Orthogonalizes the matrix.
+     * @brief Orthonormalizes the matrix.
+     *
+     * If the norm of a matrix column is below the value of zero_tol then it
+     * is considered to be zero, and we do not divide by it.
+     * Therefore, that column is considered to be zero and is not normalized.
+     * By default, zero_tol = 1.0e-15.
+     */
+    void
+    orthogonalize(double zero_tol = 1.0e-15);
+
+    /**
+     * @brief Orthonormalizes the matrix's last column, assuming the previous
+     * columns are already orthonormal.
+     *
+     * By default, ncols == -1, and the function considers the whole matrix.
+     * If ncols != -1 and ncols < d_num_cols, then a subset of the matrix
+     * is considered.
+     * This allows one to reorthonormalize the matrix every time a new column
+     * is added, assuming the previous columns have remained unchanged.
+     *
+     * If the norm of a matrix column is below the value of zero_tol then it
+     * is considered to be zero, and we do not divide by it.
+     * Therefore, that column is considered to be zero and is not normalized.
+     * By default, zero_tol = 1.0e-15.
      */
     void
-    orthogonalize();
+    orthogonalize_last(int ncols = -1, double zero_tol = 1.0e-15);
 
     /**
      * @brief Rescale every matrix row by its maximum absolute value.
diff --git a/unit_tests/test_Matrix.cpp b/unit_tests/test_Matrix.cpp
index 7534d5b19..8f3d77503 100644
--- a/unit_tests/test_Matrix.cpp
+++ b/unit_tests/test_Matrix.cpp
@@ -457,6 +457,153 @@ TEST(MatrixSerialTest, Test_Matrix_rescale_cols_max2)
                     "(i, j) = (" << i << ", " << j << ")";
 }
 
+TEST(MatrixSerialTest, Test_Matrix_orthogonalize)
+{
+    // Matrix data to orthonormalize.
+    double d_mat[16] = {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.
+    double d_mat2[16] = {1.0, 0.0, 0.0, 0.0,
+                         0.0, 1.0, 0.0, 0.0,
+                         0.0, 0.0, 1.0, 0.0,
+                         0.0, 0.0, 0.0, 1.0
+                        };
+
+    CAROM::Matrix matrix(d_mat, 4, 4, false);
+    CAROM::Matrix target(d_mat2, 4, 4, false);
+
+    double abs_error = 1.0e-15; // absolute error threshold
+
+    matrix.orthogonalize();
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 4; j++)
+            EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) <<
+                    "(i, j) = (" << i << ", " << j << ")";
+}
+
+TEST(MatrixSerialTest, Test_Matrix_orthogonalize2)
+{
+    // Matrix data to orthonormalize.
+    double d_mat[16] = {3.5, 7.1, 0.0, 0.0,
+                        0.0, 1.9, 8.3, 1e-14,
+                        0.0, 0.0, 5.7, 1.0+1.0e-14,
+                        0.0, 0.0, 0.0, 0.0
+                       };
+
+    // Target matrix data.
+    double d_mat2[16] = {1.0, 0.0, 0.0, 0.0,
+                         0.0, 1.0, 0.0, 0.0,
+                         0.0, 0.0, 1.0, 0.0,
+                         0.0, 0.0, 0.0, 0.0
+                        };
+
+    CAROM::Matrix matrix(d_mat, 4, 4, false);
+    CAROM::Matrix target(d_mat2, 4, 4, false);
+
+    double abs_error = 1.0e-15; // absolute error threshold
+
+    matrix.orthogonalize();
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 4; j++)
+            EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) <<
+                    "(i, j) = (" << i << ", " << j << ")";
+}
+
+TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last)
+{
+    // Matrix data to orthonormalize.
+    double d_mat[16] = {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.
+    double d_mat2[16] = {1.0, 0.0, 0.0, 0.0,
+                         0.0, 1.0, 0.0, 0.0,
+                         0.0, 0.0, 1.0, 0.0,
+                         0.0, 0.0, 0.0, 1.0
+                        };
+
+    CAROM::Matrix matrix(d_mat, 4, 4, false);
+    CAROM::Matrix target(d_mat2, 4, 4, false);
+
+    double abs_error = 1.0e-15; // absolute error threshold
+
+    matrix.orthogonalize_last();
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 4; j++)
+            EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) <<
+                    "(i, j) = (" << i << ", " << j << ")";
+}
+
+TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last2)
+{
+    // Matrix data to orthonormalize.
+    double d_mat[16] = {1.0, 0.0, 3.8, 1.3,
+                        0.0, 1.0, 5.6, 4.7,
+                        0.0, 0.0, 9.8, 2.5,
+                        0.0, 0.0, 0.0, 7.3
+                       };
+
+    // Target matrix data.
+    double d_mat2[16] = {1.0, 0.0, 0.0, 0.0,
+                         0.0, 1.0, 0.0, 0.0,
+                         0.0, 0.0, 1.0, 0.0,
+                         0.0, 0.0, 0.0, 1.0
+                        };
+
+    CAROM::Matrix matrix(d_mat, 4, 4, false);
+    CAROM::Matrix target(d_mat2, 4, 4, false);
+
+    double abs_error = 1.0e-15; // absolute error threshold
+
+    matrix.orthogonalize_last(3);
+    matrix.orthogonalize_last(4);
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 4; j++)
+            EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) <<
+                    "(i, j) = (" << i << ", " << j << ")";
+}
+
+TEST(MatrixSerialTest, Test_Matrix_orthogonalize_last3)
+{
+    // Matrix data to orthonormalize.
+    double d_mat[16] = {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.
+    double d_mat2[16] = {1.0, 0.0, 0.0, 0.0,
+                         0.0, 1.0, 0.0, 0.0,
+                         0.0, 0.0, 1.0, 0.0,
+                         0.0, 0.0, 0.0, 1.0
+                        };
+
+    CAROM::Matrix matrix(d_mat, 4, 4, false);
+    CAROM::Matrix target(d_mat2, 4, 4, false);
+
+    for (int i = 0; i < 4; i++)
+        matrix.orthogonalize_last(i+1);
+
+    double abs_error = 1.0e-15; // absolute error threshold
+
+    for (int i = 0; i < 4; i++)
+        for (int j = 0; j < 4; j++)
+            EXPECT_NEAR(matrix.item(i, j), target.item(i, j), abs_error) <<
+                    "(i, j) = (" << i << ", " << j << ")";
+}
+
 TEST(MatrixSerialTest, Test_pMatrix_mult_reference)
 {
     /**