generalizing to solve for multiple eigenpairs

include/lambda_lanczos/lambda_lanczos.hpp

@@ -151,7 +151,9 @@ class LambdaLanczos {
    * @brief Executes Lanczos algorithm and stores the result into reference variables passed as arguments.
    * @return Lanczos-iteration count
    */
-  int run(real_t<T>& eigvalue, std::vector<T>& eigvec) const {
+  int run(std::vector<real_t<T>>& eigvalues, std::vector<std::vector<T>>& eigvecs) const {
+    const size_t nroot = eigvecs.size();
+    assert(eigvalues.size()==nroot);
     assert(0 < this->tridiag_eps_ratio && this->tridiag_eps_ratio < 1);
 
     std::vector<std::vector<T>> u; // Lanczos vectors
@@ -178,8 +180,8 @@ class LambdaLanczos {
     util::normalize(uk);
     u.push_back(uk);
 
-    real_t<T> ev = 0.0; // Calculated eigenvalue (the initial value will never be used)
-    real_t<T> pev = std::numeric_limits<real_t<T>>::max(); // Previous eigenvalue
+    std::vector<real_t<T>> evs(nroot, 0.0); // Calculated eigenvalue (the initial value will never be used)
+    std::vector<real_t<T>> pevs(nroot, std::numeric_limits<real_t<T>>::max()); // Previous eigenvalue
 
     int itern = this->max_iteration;
     for(size_t k = 1;k <= this->max_iteration;k++) {
@@ -214,7 +216,8 @@ class LambdaLanczos {
       betak = util::norm(uk);
       beta.push_back(betak);
 
-      ev = find_mth_eigenvalue(alpha, beta, this->find_maximum ? alpha.size()-2 : 0);
+      for (size_t iroot=0ul; iroot<nroot; ++iroot)
+          evs[iroot] = find_mth_eigenvalue(alpha, beta, this->find_maximum ? alpha.size()-2-iroot : iroot);
       // The first element of alpha is a dummy. Thus its size is alpha.size()-1
 
       const real_t<T> zero_threshold = util::minimum_effective_decimal<real_t<T>>()*1e-1;
@@ -231,15 +234,24 @@ class LambdaLanczos {
       util::normalize(uk);
       u.push_back(uk);
 
-      if(std::abs(ev-pev) < std::min(std::abs(ev), std::abs(pev))*this->eps) {
-        itern = k;
-        break;
-      } else {
-        pev = ev;
+      /*
+       * only break loop if convergence condition is met for all roots
+       */
+      bool break_cond = true;
+      for (size_t iroot=0ul; iroot<nroot; ++iroot) {
+          const auto& ev = evs[iroot];
+          const auto& pev = pevs[iroot];
+          if (std::abs(ev - pev) >= std::min(std::abs(ev), std::abs(pev)) * this->eps){
+              break_cond = false;
+              break;
+          }
       }
+      if (break_cond) break;
+      else pevs = evs;
     }
 
-    eigvalue = ev - this->eigenvalue_offset;
+    eigvalues = evs;
+    for (auto& ev: eigvalues) ev -= this->eigenvalue_offset;
 
     auto m = alpha.size();
     std::vector<T> cv(m+1);
@@ -249,23 +261,27 @@ class LambdaLanczos {
 
     beta[m-1] = 0.0;
 
-    if(eigvec.size() < n) {
-      eigvec.resize(n);
-    }
+    for(size_t iroot=0ul; iroot<nroot; ++iroot) {
+        auto& eigvec = eigvecs[iroot];
+        auto& ev = evs[iroot];
 
-    for(size_t i = 0;i < n;i++) {
-      eigvec[i] = cv[m-1]*u[m-1][i];
-    }
+        if (eigvec.size() < n) {
+            eigvec.resize(n);
+        }
 
-    for(size_t k = m-2;k >= 1;k--) {
-      cv[k] = ((ev - alpha[k+1])*cv[k+1] - beta[k+1]*cv[k+2])/beta[k];
+        for (size_t i = 0; i < n; i++) {
+            eigvec[i] = cv[m - 1] * u[m - 1][i];
+        }
 
-      for(size_t i = 0;i < n;i++) {
-        eigvec[i] += cv[k]*u[k][i];
-      }
-    }
+        for (size_t k = m - 2; k >= 1; k--) {
+            cv[k] = ((ev - alpha[k + 1]) * cv[k + 1] - beta[k + 1] * cv[k + 2]) / beta[k];
 
-    util::normalize(eigvec);
+            for (size_t i = 0; i < n; i++) {
+                eigvec[i] += cv[k] * u[k][i];
+            }
+        }
+        util::normalize(eigvec);
+    }
 
     return itern;
   }
@@ -287,6 +303,15 @@ class LambdaLanczos {
     return {eigvalue, eigvec, itern};
   }
 
+  int run(real_t<T>& eigvalue, std::vector<T>& eigvec) const{
+      std::vector<real_t<T>> eigvalues(1);
+      std::vector<std::vector<T>> eigvecs(1);
+      auto retval = run(eigvalues, eigvecs);
+      eigvalue = eigvalues[0];
+      eigvec = std::move(eigvecs[0]);
+      return retval;
+  }
+
 private:
   /**
    * @brief Finds the `m`th smaller eigenvalue of given tridiagonal matrix.

src/samples/CMakeLists.txt

@@ -4,3 +4,4 @@ add_executable(sample1 sample1_simple.cpp)
 add_executable(sample2 sample2_sparse.cpp)
 add_executable(sample3 sample3_dynamic.cpp)
 add_executable(sample4 sample4_use_Eigen_library.cpp)
+add_executable(sample5 sample5_multiroot.cpp)

src/samples/sample5_multiroot.cpp

@@ -0,0 +1,57 @@
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <lambda_lanczos.hpp>
+
+using std::cout;
+using std::endl;
+using std::setprecision;
+using lambda_lanczos::LambdaLanczos;
+
+template<typename T>
+using vector = std::vector<T>;
+
+int main() {
+  const int n = 8;
+  double matrix[n][n] = {
+      { 6, -3, -3,  0, -1,  1, -1,  1},
+      {-3, -4,  2,  2, -1, -5,  0, -4},
+      {-3,  2,  2, -3,  0,  0, -1, -1},
+      { 0,  2, -3,  0, -3,  3,  2,  2},
+      {-1, -1,  0, -3, -2,  0, -5, -4},
+      { 1, -5,  0,  3,  0, -4,  5,  0},
+      {-1,  0, -1,  2, -5,  5, -4,  4},
+      { 1, -4, -1,  2, -4,  0,  4,  2}
+  };
+  /* Its eigenvalues are
+   * -13.215086, -8.500332, -4.266749, -0.467272, 0.397895, 2.303837, 7.400955, 12.346751
+   */
+
+  // the matrix-vector multiplication routine
+  auto mv_mul = [&](const vector<double>& in, vector<double>& out) {
+    for(int i = 0;i < n;i++) {
+      for(int j = 0;j < n;j++) {
+        out[i] += matrix[i][j]*in[j];
+      }
+    } 
+  };
+
+  LambdaLanczos<double> engine(mv_mul, n, false); // true means to calculate the largest eigenvalue.
+
+  const size_t nroot = 2;
+  vector<double> eigenvalues(nroot);
+  vector<vector<double>> eigenvectors(nroot, {n});
+  int itern = engine.run(eigenvalues, eigenvectors);
+
+    cout << "Iteration count: " << itern << endl;
+  for (size_t iroot=0ul; iroot<nroot; ++iroot) {
+      cout << "Eigen value (root "<< iroot <<"): " << setprecision(16) << eigenvalues[iroot] << endl;
+      cout << "Eigen vector (root "<< iroot <<"): ";
+      for (int i = 0; i < n; i++) {
+          cout << eigenvectors[iroot][i] << " ";
+      }
+      cout << endl;
+  }
+
+  return EXIT_SUCCESS;
+}