allows non square matrix for the mgs

src/algorithms/davidson.rs

@@ -129,7 +129,6 @@ impl EigenDavidson {
                 result = Ok(create_results(&eig.eigenvalues, &ritz_vectors, nvalues));
                 break;
             }
-
             // 5. Update subspace basis set
             // 5.1 Add the correction vectors to the current basis
             if dim_sub + nvalues <= conf.max_search_space {
@@ -138,7 +137,8 @@ impl EigenDavidson {
                 update_subspace(&mut basis, correction, dim_sub, dim_sub + nvalues);
 
                 // 6. Orthogonalize the subspace
-                basis = orthogonalize_subspace(basis);
+                basis = orthogonalize_subspace(basis, dim_sub, dim_sub + nvalues);
+
                 // update counter
                 dim_sub += nvalues;
 
@@ -269,8 +269,8 @@ fn update_subspace(basis: &mut DMatrix<f64>, vectors: DMatrix<f64>, start: usize
 }
 
 /// Orthogonalize the subpsace using the QR method
-fn orthogonalize_subspace(vectors: DMatrix<f64>) -> DMatrix<f64> {
-    let mgs = MGS::new(vectors);
+fn orthogonalize_subspace(vectors: DMatrix<f64>, start: usize, end: usize) -> DMatrix<f64> {
+    let mgs = MGS::new(vectors, start, end);
     match mgs {
         Ok(result) => result.basis,
         Err(msg) => panic!("Error orthonormalising the basis:{}", msg),

src/modified_gram_schmidt.rs

@@ -16,19 +16,13 @@ pub struct MGS {
 }
 
 impl MGS {
-    /// The new static method takes a single argument:
+    /// The new static method takes three argument:
     /// * `vectors` to diagonalize as columns of the matrix
-    pub fn new(vectors: DMatrix<f64>) -> Result<Self, &'static str> {
-        if !vectors.is_square() {
-            return Err("vectors must be an square matrix");
-        }
-        let dim = vectors.nrows();
-        let mut basis = DMatrix::<f64>::zeros(dim, dim);
-        // first vector of the basis
-        let first = vectors.column(0) / vectors.column(0).norm();
-        basis.set_column(0, &first);
-        // Iterate over the rest of the columns
-        for i in 1..basis.ncols() {
+    /// * `start` index of the column to start orthogonalizing
+    /// * `end` last index of the column to diagonalize (non-inclusive)
+    pub fn new(vectors: DMatrix<f64>, start: usize, end: usize) -> Result<Self, &'static str> {
+        let mut basis = vectors.clone();
+        for i in start..end {
             basis.set_column(i, &vectors.column(i));
             for j in 0..i {
                 let proj = MGS::project(&basis.column(j), &basis.column(i));
@@ -55,12 +49,22 @@ mod test {
     fn test_gram_schmidt() {
         let dim = 10;
         let vectors = DMatrix::<f64>::new_random(dim, dim);
-        let mgs_result = super::MGS::new(vectors);
-        let basis: DMatrix<f64> = match mgs_result {
+        fun_test(vectors, 0);
+    }
+
+    #[test]
+    fn test_start_gram_schmidt() {
+        let arr = DMatrix::<f64>::from_vec(3, 3, vec![1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 2.0]);
+        fun_test(arr, 1);
+    }
+
+    fn fun_test(vectors: DMatrix<f64>, start: usize) {
+        let end = vectors.ncols();
+        let mgs_result = super::MGS::new(vectors, start, end);
+        let basis = match mgs_result {
             Ok(ortho) => ortho.basis,
             Err(message) => panic!(message),
         };
-
         let result = basis.transpose() * &basis;
         assert!(result.is_identity(1e-8));
     }

tests/tests_davidson.rs

@@ -35,15 +35,17 @@ fn test_davidson_unsorted() {
 #[test]
 fn test_davidson_highest() {
     // Test the compution of the highest eigenvalues
-    let arr = generate_diagonal_dominant(10, 0.005);
+    let dim = 100;
+    let nvalues = 7;
+    let arr = generate_diagonal_dominant(dim, 0.005);
     let eig = sort_eigenpairs(na::linalg::SymmetricEigen::new(arr.clone()), false);
 
     let target = SpectrumTarget::Highest;
 
-    let dav_eig = EigenDavidson::new(arr.clone(), 2, "DPR", target.clone()).unwrap();
-    test_eigenpairs(&eig, dav_eig, 2);
-    let dav_eig = EigenDavidson::new(arr.clone(), 2, "GJD", target).unwrap();
-    test_eigenpairs(&eig, dav_eig, 2);
+    let dav_eig = EigenDavidson::new(arr.clone(), nvalues, "DPR", target.clone()).unwrap();
+    test_eigenpairs(&eig, dav_eig, nvalues);
+    let dav_eig = EigenDavidson::new(arr.clone(), nvalues, "GJD", target).unwrap();
+    test_eigenpairs(&eig, dav_eig, nvalues);
 }
 
 fn test_eigenpairs(