Implement DPR

include/Spectra/JDSymEigsBase.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2016-2020 Yixuan Qiu <[email protected]>
+// Copyright (C) 2020 Netherlands eScience Center <[email protected]>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -38,7 +38,6 @@ class JDSymEigsBase
 public:
     JDSymEigsBase(OpType& op, Index nev) :
         matrix_operator_(op),
-        operator_dimension_(op.rows()),
         number_eigenvalues_(nev),
         max_search_space_size_(10 * number_eigenvalues_),
         initial_search_space_size_(2 * number_eigenvalues_),
@@ -92,27 +91,27 @@ class JDSymEigsBase
     Matrix eigenvectors() const { return ritz_pairs_.RitzVectors().leftCols(number_eigenvalues_); }
 
 protected:
-    virtual Matrix SetupInitialSearchSpace() const = 0;
+    virtual Matrix SetupInitialSearchSpace(SortRule selection) const = 0;
 
     virtual Matrix CalculateCorrectionVector() const = 0;
-    const OpType& matrix_operator_;  // object to conduct marix operation,
+    const OpType& matrix_operator_;  // object to conduct matrix operation,
                                      // e.g. matrix-vector product
 
     Index niter_ = 0;
-    const Index operator_dimension_;  // dimension of matrix A
     const Index number_eigenvalues_;  // number of eigenvalues requested
     Index max_search_space_size_;
     Index initial_search_space_size_;
     Index correction_size_;             // how many correction vectors are added in each iteration
     RitzPairs<Scalar> ritz_pairs_;      // Ritz eigen pair structure
     SearchSpace<Scalar> search_space_;  // search space
+    Index size_update_;                 // size of the current correction
 
 private:
     CompInfo info_ = CompInfo::NotComputed;  // status of the computation
 
     void check_argument() const
     {
-        if (number_eigenvalues_ < 1 || number_eigenvalues_ > operator_dimension_ - 1)
+        if (number_eigenvalues_ < 1 || number_eigenvalues_ > matrix_operator_.cols() - 1)
             throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 1, n is the size of matrix");
     }
 

include/Spectra/JDSymEigsDPR.h

@@ -0,0 +1,88 @@
+// Copyright (C) 2020 Netherlands eScience Center <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#ifndef SPECTRA_JD_SYM_EIGS_DPR_H
+#define SPECTRA_JD_SYM_EIGS_DPR_H
+
+#include <Eigen/Dense>
+#include "JDSymEigsBase.h"
+#include "Util/SelectionRule.h"
+
+namespace Spectra {
+
+///
+/// \ingroup EigenSolver
+///
+/// This class implement the DPR correction for the Davidson algorithms.
+/// The algorithms in the Davidson family only differ in how the correction
+/// vectors are computed and optionally in the initial orthogonal basis set.
+///
+/// the DPR correction compute the new correction vector using the following expression:
+/// \f[ correction = -(\boldsymbol{D} - \rho \boldsymbol{I})^{-1} \boldsymbol{r} \f]
+/// where
+/// \f$D\f$ is the diagonal of the target matrix, \f$\rho\f$ the Ritz eigenvalue,
+/// \f$I\f$ the identity matrix and \f$r\f$ the residue vector.
+///
+
+template <typename OpType>
+class JDSymEigsDPR : public JDSymEigsBase<OpType>
+{
+private:
+    using Index = Eigen::Index;
+    using Scalar = typename OpType::Scalar;
+    using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+    Vector diagonal_;
+
+public:
+    JDSymEigsDPR(OpType& op, Index nev) :
+        JDSymEigsBase<OpType>{op, nev}
+    {
+        diagonal_.resize(this->matrix_operator_.rows());
+        for (Index i = 0; i < op.rows(); i++)
+        {
+            diagonal_(i) = op(i, i);
+        }
+    }
+
+    /// Create initial search space based on the diagonal
+    /// and the spectrum'target (highest or lowest)
+    /// \param selection spectrum section to target (e.g. lowest, etc.)
+    /// \return Matrix with the initial orthonormal basis
+    Matrix SetupInitialSearchSpace(SortRule selection) const final
+    {
+        std::vector<Eigen::Index> indices_sorted = argsort(selection, diagonal_);
+
+        Matrix initial_basis = Matrix::Zero(this->matrix_operator_.rows(), this->initial_search_space_size_);
+
+        for (Index k = 0; k < this->initial_search_space_size_; k++)
+        {
+            Index row = indices_sorted[k];
+            initial_basis(row, k) = 1.0;
+        }
+        return initial_basis;
+    }
+
+    /// compute the corrections using the DPR method.
+    /// \return new correction vectors.
+    Matrix CalculateCorrectionVector() const final
+    {
+        const Matrix& residues = this->ritz_pairs_.Residues();
+        const Vector& eigvals = this->ritz_pairs_.RitzValues();
+        Matrix correction = Matrix::Zero(this->matrix_operator_.rows(), this->correction_size_);
+        for (Index k = 0; k < this->correction_size_; k++)
+        {
+            Vector tmp = eigvals(k) - diagonal_.array();
+            correction.col(k) = residues.col(k).array() / tmp.array();
+        }
+        return correction;
+    }
+};
+
+}  // namespace Spectra
+
+#endif  // SPECTRA_JD_SYM_EIGS_DPR_H

include/Spectra/Util/SearchSpace.h

@@ -25,7 +25,7 @@ class SearchSpace
 public:
     SearchSpace() = default;
 
-    void size() const
+    Index size() const
     {
         return basis_vectors_.cols();
     }

test/CMakeLists.txt

@@ -8,18 +8,18 @@ list(APPEND test_target_sources
         GenEigs.cpp
         GenEigsRealShift.cpp
         GenEigsComplexShift.cpp
+        JDSymEigsBase.cpp
+        JDSymEigsDPR.cpp
         QR.cpp
+        RitzPairs.cpp
+        SearchSpace.cpp
         SparseGenMatProd.cpp
         SparseSymMatProd.cpp
         SVD.cpp
         SymEigs.cpp
         SymEigsShift.cpp
         SymGEigsCholesky.cpp
         SymGEigsRegInv.cpp
-        SymGEigsShift.cpp
-        JDSymEigsBase.cpp
-        RitzPairs.cpp
-        SearchSpace.cpp
         )
 
 foreach(TEST_SOURCE ${test_target_sources})

test/JDSymEigsBase.cpp

@@ -23,7 +23,7 @@ class JDMock : public JDSymEigsBase<OpType>
 public:
     JDMock(OpType& op, Index nev) :
         JDSymEigsBase<OpType>(op, nev) {}
-    Matrix SetupInitialSearchSpace() const
+    Matrix SetupInitialSearchSpace(SortRule selection) const
     {
         return Matrix::Zero(0, 0);
     }

test/JDSymEigsDPR.cpp

@@ -0,0 +1,20 @@
+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+
+#include <Spectra/JDSymEigsDPR.h>
+#include <Spectra/MatOp/DenseGenMatProd.h>
+
+using namespace Spectra;
+
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+TEMPLATE_TEST_CASE("Constructing JDSymEigsDPR", "[JDSymEigsDPR]", float, double)
+{
+    using Matrix = Eigen::Matrix<TestType, Eigen::Dynamic, Eigen::Dynamic>;
+    const Matrix A = Matrix::Random(10, 10);
+    DenseGenMatProd<TestType> op(A);
+    JDSymEigsDPR<DenseGenMatProd<TestType>> eigs{op, 5};
+    REQUIRE(eigs.num_iterations() == 0);
+    REQUIRE(eigs.info() == Spectra::CompInfo::NotComputed);
+}