fixed DPR implementation #21

include/Spectra/JDSymEigsDPR.h

@@ -17,14 +17,18 @@ template <typename Scalar, typename OpType>
 class JDSymEigsDPR : public JDSymEigsBase<Scalar, OpType>
 {
 private:
-    Vector diagonal_(operator_dimension_);
+    using Index = Eigen::Index;
+    using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+    Vector diagonal_ = Vector::Zero(this->operator_dimension_);
     std::vector<Eigen::Index> indices_sorted_;
 
     void extract_diagonal()
     {
-        for (Index i = 0; i < operator_dimension_; i++)
+        for (Index i = 0; i < this->operator_dimension_; i++)
         {
-            diagonal_(i) = matrix_operator_(i, i);
+            diagonal_(i) = this->matrix_operator_(i, i);
         }
     }
 
@@ -41,12 +45,12 @@ class JDSymEigsDPR : public JDSymEigsBase<Scalar, OpType>
         extract_diagonal();
         calculate_indices_diagonal_sorted(selection);
 
-        Matrix initial_basis = Matrix::Zero(operator_dimension_, initial_search_space_size_);
+        Matrix initial_basis = Matrix::Zero(this->operator_dimension_, this->initial_search_space_size_);
 
-        for (Index k = 0; k < initial_search_space_size_; k++)
+        for (Index k = 0; k < this->initial_search_space_size_; k++)
         {
             Index row = indices_sorted_[k];
-            initial_basis.coeff(row, k) = diagonal_[k];
+            initial_basis.coeff(row, k) = 1.0;
         }
         return initial_basis;
     }
@@ -55,13 +59,12 @@ class JDSymEigsDPR : public JDSymEigsBase<Scalar, OpType>
     /// \return new correction vectors.
     Matrix CalculateCorrectionVector() const final
     {
-        Index nresidues = search_space_.size();
-        const Matrix& residues = Residues();
-        const Vector& eigenvalues = ritz_pairs_.RitzValues();
-        Matrix correction = Matrix::zero(operator_dimension_, nresidues);
-        for (Index k = 0; k < nresidues; k++)
+        const Matrix& residues = this->Residues();
+        const Vector& eigvals = this->ritz_pairs_.RitzValues();
+        Matrix correction = Matrix::Zero(this->operator_dimension_, this->size_update_);
+        for (Index k = 0; k < this->size_update_; k++)
         {
-            Vector tmp = Vector::Constant(eigenvalues(k)) - diagonal_;
+            Vector tmp = eigvals(k) - diagonal_.array();
             correction.col(k) = residues.col(k).array() / tmp.array();
         }
         return correction;

test/JDSymEigsDPR.cpp

@@ -0,0 +1,14 @@
+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+
+#include <Spectra/JDSymEigsDPR.h>
+
+using namespace Spectra;
+
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+TEST_CASE("Constructing JDSymEigsDPR", "[JDSymEigsDPR]")
+{
+    REQUIRE(true);
+}