From 4c379205547bdc5bfa945bede31db6d52c113340 Mon Sep 17 00:00:00 2001
From: felipez <tifonzafel@gmail.com>
Date: Wed, 1 Jul 2020 12:55:40 +0200
Subject: [PATCH] fixed DPR implementation #21

---
 include/Spectra/JDSymEigsDPR.h | 27 +++++++++++++++------------
 test/JDSymEigsDPR.cpp          | 14 ++++++++++++++
 2 files changed, 29 insertions(+), 12 deletions(-)
 create mode 100644 test/JDSymEigsDPR.cpp

diff --git a/include/Spectra/JDSymEigsDPR.h b/include/Spectra/JDSymEigsDPR.h
index f4a24e5e..50e90fee 100644
--- a/include/Spectra/JDSymEigsDPR.h
+++ b/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;
diff --git a/test/JDSymEigsDPR.cpp b/test/JDSymEigsDPR.cpp
new file mode 100644
index 00000000..76f75511
--- /dev/null
+++ b/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);
+}
\ No newline at end of file