used vectorize modified GS #11

NLESC-JCER · Jun 30, 2020 · 539cfc3 · 539cfc3
1 parent 86cd4d1
commit 539cfc3
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 11 deletions.
diff --git a/include/Spectra/LinAlg/Orthogonalization.h b/include/Spectra/LinAlg/Orthogonalization.h
@@ -54,11 +54,11 @@ class Orthogonalization
     Orthogonalization(const Matrix& A) :
         basis{A} {}
 
-    /// two consecutive Gram-schmidt iterations are enough to converge
+    /// two consecutive modified Gram-schmidt iterations are enough to converge
     /// http://stoppels.blog/posts/orthogonalization-performance
     /// \return Returned matrix type will be `Eigen::Matrix<Scalar, ...>`, depending on
     /// the template parameter `Scalar` defined.
-    Matrix twice_gramschmidt()
+    Matrix twice_modified_gramschmidt()
     {
         Matrix Q = basis;
         Index nstart = start_index;
@@ -92,12 +92,7 @@ class Orthogonalization
         }
         for (Index j = nstart; j < basis.cols(); ++j)
         {
-            Vector tmp = Q.col(j);
-            for (Index k = 0; k < j; k++)
-            {
-                tmp -= Q.col(k).dot(tmp) * Q.col(k);
-            }
-            Q.col(j) = tmp;
+            Q.col(j) -= Q.leftCols(j) * (Q.leftCols(j).transpose() * Q.col(j));
             check_linear_dependency(Q.col(j), j);
             Q.col(j).normalize();
         }

diff --git a/test/Orthogonalization.cpp b/test/Orthogonalization.cpp
@@ -26,7 +26,7 @@ TEST_CASE("complete orthonormalization", "[Gram-Schmidt]")
 
     MatrixXd mat = MatrixXd::Random(n, n);
     Orthogonalization<double> gs{mat};
-    check_orthogonality(gs.twice_gramschmidt());
+    check_orthogonality(gs.twice_modified_gramschmidt());
     check_orthogonality(gs.modified_gramschmidt());
     check_orthogonality(gs.QR());
 }
@@ -39,14 +39,14 @@ TEST_CASE("Partial orthonormalization", "[Gram-Schmidt]")
     // Create a n x 20 orthonormal basis
     MatrixXd mat = MatrixXd::Random(n, n - 20);
     Orthogonalization<double> gs{mat};
-    mat.leftCols(n - 20) = gs.twice_gramschmidt();
+    mat.leftCols(n - 20) = gs.twice_modified_gramschmidt();
 
     mat.conservativeResize(Eigen::NoChange, n);
     mat.rightCols(20) = MatrixXd::Random(n, 20);
 
     // Orthogonalize from 80 onwards
     Orthogonalization<double> new_gs{mat, 80};
-    check_orthogonality(new_gs.twice_gramschmidt());
+    check_orthogonality(new_gs.twice_modified_gramschmidt());
     check_orthogonality(new_gs.modified_gramschmidt());
     check_orthogonality(new_gs.QR());
 }