Logging

include/Spectra/DavidsonSymEigsSolver.h

@@ -52,6 +52,19 @@ class DavidsonSymEigsSolver : public JDSymEigsBase<DavidsonSymEigsSolver<OpType>
     DavidsonSymEigsSolver(OpType& op, Index nev) :
         DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev) {}
 
+    DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, LoggerBase<Scalar, Vector>* logger) :
+        JDSymEigsBase<DavidsonSymEigsSolver<OpType>, OpType>(op, nev, nvec_init, nvec_max, logger)
+    {
+        m_diagonal.resize(this->m_matrix_operator.rows());
+        for (Index i = 0; i < op.rows(); i++)
+        {
+            m_diagonal(i) = op(i, i);
+        }
+    }
+
+    DavidsonSymEigsSolver(OpType& op, Index nev, LoggerBase<Scalar, Vector>* logger) :
+        DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, logger) {}
+
     /// Create initial search space based on the diagonal
     /// and the spectrum'target (highest or lowest)
     ///

include/Spectra/GenEigsBase.h

@@ -24,6 +24,7 @@
 #include "LinAlg/DoubleShiftQR.h"
 #include "LinAlg/UpperHessenbergEigen.h"
 #include "LinAlg/Arnoldi.h"
+#include "LoggerBase.h"
 
 namespace Spectra {
 
@@ -70,10 +71,12 @@ class GenEigsBase
     ComplexVector m_ritz_val;  // Ritz values
     ComplexMatrix m_ritz_vec;  // Ritz vectors
     ComplexVector m_ritz_est;  // last row of m_ritz_vec, also called the Ritz estimates
+    ComplexVector m_resid;
 
 private:
     BoolArray     m_ritz_conv; // indicator of the convergence of Ritz values
     CompInfo      m_info;      // status of the computation
+    LoggerBase<Scalar, ComplexVector> *m_logger = nullptr;
     // clang-format on
 
     // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part
@@ -150,9 +153,9 @@ class GenEigsBase
 
         // thresh = tol * max(eps23, abs(theta)), theta for Ritz value
         Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(eps23);
-        Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac.f_norm();
+        m_resid = m_ritz_est.head(m_nev).array().abs().matrix() * m_fac.f_norm();
         // Converged "wanted" Ritz values
-        m_ritz_conv = (resid < thresh);
+        m_ritz_conv = (m_resid.real().array() < thresh);
 
         return m_ritz_conv.count();
     }
@@ -321,7 +324,6 @@ class GenEigsBase
 
 public:
     /// \cond
-
     GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv) :
         m_op(op),
         m_n(m_op.rows()),
@@ -339,6 +341,19 @@ class GenEigsBase
             throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix");
     }
 
+    GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, LoggerBase<Scalar, ComplexVector>* logger) :
+        GenEigsBase(op, Bop, nev, ncv)
+    {
+        set_logger(logger);
+    }
+
+    ///
+    /// Sets the logger ptr with a user constructed logger object.
+    ///
+    void set_logger(LoggerBase<Scalar, ComplexVector>* logger)
+    {
+        m_logger = logger;
+    }
     ///
     /// Virtual destructor
     ///
@@ -362,11 +377,13 @@ class GenEigsBase
         m_ritz_vec.resize(m_ncv, m_nev);
         m_ritz_est.resize(m_ncv);
         m_ritz_conv.resize(m_nev);
+        m_resid.resize(m_nev);
 
         m_ritz_val.setZero();
         m_ritz_vec.setZero();
         m_ritz_est.setZero();
         m_ritz_conv.setZero();
+        m_resid.setZero();
 
         m_nmatop = 0;
         m_niter = 0;
@@ -424,12 +441,29 @@ class GenEigsBase
         Index i, nconv = 0, nev_adj;
         for (i = 0; i < maxit; i++)
         {
+            if (m_logger)
+            {
+                m_logger->call_iteration_start();
+            }
             nconv = num_converged(tol);
+            const ComplexVector eigs = m_ritz_val.head(m_nev);
+            const IterationData<Scalar, ComplexVector> data(i, nconv, m_ncv, eigs, m_resid, m_ritz_conv);
+
             if (nconv >= m_nev)
+            {
+                if (m_logger)
+                {
+                    m_logger->call_iteration_end(data);
+                }
                 break;
+            }
 
             nev_adj = nev_adjusted(nconv);
             restart(nev_adj, selection);
+            if (m_logger)
+            {
+                m_logger->call_iteration_end(data);
+            }
         }
         // Sorting results
         sort_ritzpair(sorting);

include/Spectra/GenEigsComplexShiftSolver.h

@@ -37,6 +37,7 @@ class GenEigsComplexShiftSolver : public GenEigsBase<OpType, IdentityBOp>
     using Index = Eigen::Index;
     using Complex = std::complex<Scalar>;
     using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+    using ComplexVector = Eigen::Matrix<Complex, Eigen::Dynamic, 1>;
     using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
 
     using Base = GenEigsBase<OpType, IdentityBOp>;
@@ -152,6 +153,32 @@ class GenEigsComplexShiftSolver : public GenEigsBase<OpType, IdentityBOp>
     {
         op.set_shift(m_sigmar, m_sigmai);
     }
+    ///
+    /// Constructor to create a eigen solver object using the shift-and-invert mode with logging.
+    ///
+    /// \param op      The matrix operation object that implements
+    ///                the complex shift-solve operation of \f$A\f$: calculating
+    ///                \f$\mathrm{Re}\{(A-\sigma I)^{-1}v\}\f$ for any vector \f$v\f$. Users could either
+    ///                create the object from the wrapper class such as DenseGenComplexShiftSolve, or
+    ///                define their own that implements all the public members
+    ///                as in DenseGenComplexShiftSolve.
+    /// \param nev     Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
+    ///                where \f$n\f$ is the size of matrix.
+    /// \param ncv     Parameter that controls the convergence speed of the algorithm.
+    ///                Typically a larger `ncv` means faster convergence, but it may
+    ///                also result in greater memory use and more matrix operations
+    ///                in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
+    ///                and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
+    /// \param sigmar  The real part of the shift.
+    /// \param sigmai  The imaginary part of the shift.
+    /// \param logger  A logging object that inherits from the base class in LoggerBase.h
+    ///
+    GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai, LoggerBase<Scalar, ComplexVector>* logger) :
+        Base(op, IdentityBOp(), nev, ncv, logger),
+        m_sigmar(sigmar), m_sigmai(sigmai)
+    {
+        op.set_shift(m_sigmar, m_sigmai);
+    }
 };
 
 }  // namespace Spectra

include/Spectra/GenEigsRealShiftSolver.h

@@ -36,6 +36,7 @@ class GenEigsRealShiftSolver : public GenEigsBase<OpType, IdentityBOp>
     using Scalar = typename OpType::Scalar;
     using Index = Eigen::Index;
     using Complex = std::complex<Scalar>;
+    using ComplexVector = Eigen::Matrix<Complex, Eigen::Dynamic, 1>;
     using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
 
     using Base = GenEigsBase<OpType, IdentityBOp>;
@@ -78,6 +79,31 @@ class GenEigsRealShiftSolver : public GenEigsBase<OpType, IdentityBOp>
     {
         op.set_shift(m_sigma);
     }
+    ///
+    /// Constructor to create a eigen solver object using the shift-and-invert mode with logging.
+    ///
+    /// \param op     The matrix operation object that implements
+    ///               the shift-solve operation of \f$A\f$: calculating
+    ///               \f$(A-\sigma I)^{-1}v\f$ for any vector \f$v\f$. Users could either
+    ///               create the object from the wrapper class such as DenseGenRealShiftSolve, or
+    ///               define their own that implements all the public members
+    ///               as in DenseGenRealShiftSolve.
+    /// \param nev    Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
+    ///               where \f$n\f$ is the size of matrix.
+    /// \param ncv    Parameter that controls the convergence speed of the algorithm.
+    ///               Typically a larger `ncv` means faster convergence, but it may
+    ///               also result in greater memory use and more matrix operations
+    ///               in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
+    ///               and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
+    /// \param sigma  The real-valued shift.
+    /// \param logger A logging object that inherits from the base class in LoggerBase.h
+    ///
+    GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, LoggerBase<Scalar, ComplexVector>* logger) :
+        Base(op, IdentityBOp(), nev, ncv, logger),
+        m_sigma(sigma)
+    {
+        op.set_shift(m_sigma);
+    }
 };
 
 }  // namespace Spectra

include/Spectra/GenEigsSolver.h

@@ -120,6 +120,9 @@ class GenEigsSolver : public GenEigsBase<OpType, IdentityBOp>
 {
 private:
     using Index = Eigen::Index;
+    using Scalar = typename OpType::Scalar;
+    using Complex = std::complex<Scalar>;
+    using ComplexVector = Eigen::Matrix<Complex, Eigen::Dynamic, 1>;
 
 public:
     ///
@@ -142,6 +145,27 @@ class GenEigsSolver : public GenEigsBase<OpType, IdentityBOp>
     GenEigsSolver(OpType& op, Index nev, Index ncv) :
         GenEigsBase<OpType, IdentityBOp>(op, IdentityBOp(), nev, ncv)
     {}
+    ///
+    /// Constructor to create a solver object with logging.
+    ///
+    /// \param op   The matrix operation object that implements
+    ///             the matrix-vector multiplication operation of \f$A\f$:
+    ///             calculating \f$Av\f$ for any vector \f$v\f$. Users could either
+    ///             create the object from the wrapper class such as DenseGenMatProd, or
+    ///             define their own that implements all the public members
+    ///             as in DenseGenMatProd.
+    /// \param nev  Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
+    ///             where \f$n\f$ is the size of matrix.
+    /// \param ncv  Parameter that controls the convergence speed of the algorithm.
+    ///             Typically a larger `ncv` means faster convergence, but it may
+    ///             also result in greater memory use and more matrix operations
+    ///             in each iteration. This parameter must satisfy \f$nev+2 \le ncv \le n\f$,
+    ///             and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
+    /// \param logger A logging object that inherits from the base class in LoggerBase.h
+    ///
+    GenEigsSolver(OpType& op, Index nev, Index ncv, LoggerBase<Scalar, ComplexVector>* logger) :
+        GenEigsBase<OpType, IdentityBOp>(op, IdentityBOp(), nev, ncv, logger)
+    {}
 };
 
 }  // namespace Spectra

include/Spectra/JDSymEigsBase.h

@@ -12,12 +12,14 @@
 #include <cmath>      // std::abs, std::pow
 #include <algorithm>  // std::min
 #include <stdexcept>  // std::invalid_argument
+#include <utility>    // std::move
 #include <iostream>
 
 #include "Util/SelectionRule.h"
 #include "Util/CompInfo.h"
 #include "LinAlg/SearchSpace.h"
 #include "LinAlg/RitzPairs.h"
+#include "LoggerBase.h"
 
 namespace Spectra {
 
@@ -51,6 +53,7 @@ class JDSymEigsBase
 
 private:
     CompInfo m_info = CompInfo::NotComputed;  // status of the computation
+    LoggerBase<Scalar, Vector>* m_logger = nullptr;
 
     void check_argument() const
     {
@@ -84,6 +87,26 @@ class JDSymEigsBase
         initialize();
     }
 
+    JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, LoggerBase<Scalar, Vector>* logger) :
+        m_matrix_operator(op),
+        m_number_eigenvalues(nev),
+        m_max_search_space_size(nvec_max < op.rows() ? nvec_max : 10 * m_number_eigenvalues),
+        m_initial_search_space_size(nvec_init < op.rows() ? nvec_init : 2 * m_number_eigenvalues),
+        m_correction_size(m_number_eigenvalues)
+    {
+        check_argument();
+        initialize();
+        set_logger(logger);
+    }
+
+    ///
+    /// Sets the logger ptr with a user constructed logger object.
+    ///
+    void set_logger(LoggerBase<Scalar, Vector>* logger)
+    {
+        m_logger = logger;
+    }
+
     JDSymEigsBase(OpType& op, Index nev) :
         JDSymEigsBase(op, nev, 2 * nev, 10 * nev) {}
 
@@ -148,6 +171,11 @@ class JDSymEigsBase
         niter_ = 0;
         for (niter_ = 0; niter_ < maxit; niter_++)
         {
+            if (m_logger)
+            {
+                m_logger->call_iteration_start();
+            }
+
             bool do_restart = (m_search_space.size() > m_max_search_space_size);
 
             if (do_restart)
@@ -166,20 +194,39 @@ class JDSymEigsBase
             m_ritz_pairs.sort(selection);
 
             bool converged = m_ritz_pairs.check_convergence(tol, m_number_eigenvalues);
+            const Eigen::Array<bool, Eigen::Dynamic, 1> conv_eig = m_ritz_pairs.converged_eigenvalues().head(m_number_eigenvalues);
+            const Index num_conv = conv_eig.count();
+            const Vector evals = eigenvalues();
+            const Vector res = m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues);
+            const Index search_space_size = m_search_space.size();
+            const IterationData<Scalar, Vector> data(niter_, num_conv, search_space_size, evals, res, conv_eig);
+
             if (converged)
             {
                 m_info = CompInfo::Successful;
+                if (m_logger)
+                {
+                    m_logger->call_iteration_end(data);
+                }
                 break;
             }
             else if (niter_ == maxit - 1)
             {
                 m_info = CompInfo::NotConverging;
+                if (m_logger)
+                {
+                    m_logger->call_iteration_end(data);
+                }
                 break;
             }
             Derived& derived = static_cast<Derived&>(*this);
             Matrix corr_vect = derived.calculate_correction_vector();
 
             m_search_space.extend_basis(corr_vect);
+            if (m_logger)
+            {
+                m_logger->call_iteration_end(data);
+            }
         }
         return (m_ritz_pairs.converged_eigenvalues()).template cast<Index>().head(m_number_eigenvalues).sum();
     }