From 32bd4939ae670c7168a84eed973f74d89dcc39b9 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Thu, 19 Aug 2021 10:23:40 -0400 Subject: [PATCH 01/16] GenEigs* logging --- include/Spectra/GenEigsBase.h | 24 +++- include/Spectra/GenEigsComplexShiftSolver.h | 8 +- include/Spectra/GenEigsRealShiftSolver.h | 8 +- include/Spectra/GenEigsSolver.h | 10 +- include/Spectra/LoggerBase.h | 47 ++++++++ test/CMakeLists.txt | 1 + test/LoggingGenEigs.cpp | 127 ++++++++++++++++++++ 7 files changed, 213 insertions(+), 12 deletions(-) create mode 100644 include/Spectra/LoggerBase.h create mode 100644 test/LoggingGenEigs.cpp diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 2a512bb..413c8e6 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -13,6 +13,7 @@ #include // std::min, std::copy #include // std::complex, std::conj, std::norm, std::abs #include // std::invalid_argument +#include // std::unique_ptr #include "Util/Version.h" #include "Util/TypeTraits.h" @@ -24,6 +25,7 @@ #include "LinAlg/DoubleShiftQR.h" #include "LinAlg/UpperHessenbergEigen.h" #include "LinAlg/Arnoldi.h" +#include "LoggerBase.h" namespace Spectra { @@ -70,10 +72,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 + std::unique_ptr> m_logger; // clang-format on // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part @@ -150,9 +154,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(); } @@ -322,7 +326,7 @@ class GenEigsBase public: /// \cond - GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv) : + GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : m_op(op), m_n(m_op.rows()), m_nev(nev), @@ -337,8 +341,18 @@ class GenEigsBase if (ncv < nev + 2 || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix"); + + set_logger(logger); } + /// + /// Sets the logger unique_ptr with a user constructed logger object. + /// + void set_logger(std::unique_ptr>& logger) + { + if (logger) + m_logger = std::move(logger); + } /// /// Virtual destructor /// @@ -362,11 +376,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; @@ -425,6 +441,8 @@ class GenEigsBase for (i = 0; i < maxit; i++) { nconv = num_converged(tol); + if (m_logger) + m_logger->iteration_log(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); if (nconv >= m_nev) break; diff --git a/include/Spectra/GenEigsComplexShiftSolver.h b/include/Spectra/GenEigsComplexShiftSolver.h index 098f6d5..1b1e692 100644 --- a/include/Spectra/GenEigsComplexShiftSolver.h +++ b/include/Spectra/GenEigsComplexShiftSolver.h @@ -37,6 +37,7 @@ class GenEigsComplexShiftSolver : public GenEigsBase using Index = Eigen::Index; using Complex = std::complex; using Vector = Eigen::Matrix; + using ComplexVector = Eigen::Matrix; using ComplexArray = Eigen::Array; using Base = GenEigsBase; @@ -128,7 +129,7 @@ class GenEigsComplexShiftSolver : public GenEigsBase public: /// - /// Constructor to create a eigen solver object using the shift-and-invert mode. + /// 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 @@ -145,9 +146,10 @@ class GenEigsComplexShiftSolver : public GenEigsBase /// 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) : - Base(op, IdentityBOp(), nev, ncv), + GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai, std::unique_ptr> logger = nullptr) : + Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigmar(sigmar), m_sigmai(sigmai) { op.set_shift(m_sigmar, m_sigmai); diff --git a/include/Spectra/GenEigsRealShiftSolver.h b/include/Spectra/GenEigsRealShiftSolver.h index 6d8dfcb..52a19af 100644 --- a/include/Spectra/GenEigsRealShiftSolver.h +++ b/include/Spectra/GenEigsRealShiftSolver.h @@ -36,6 +36,7 @@ class GenEigsRealShiftSolver : public GenEigsBase using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Complex = std::complex; + using ComplexVector = Eigen::Matrix; using ComplexArray = Eigen::Array; using Base = GenEigsBase; @@ -55,7 +56,7 @@ class GenEigsRealShiftSolver : public GenEigsBase public: /// - /// Constructor to create a eigen solver object using the shift-and-invert mode. + /// 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 @@ -71,9 +72,10 @@ class GenEigsRealShiftSolver : public GenEigsBase /// 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) : - Base(op, IdentityBOp(), nev, ncv), + GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigma(sigma) { op.set_shift(m_sigma); diff --git a/include/Spectra/GenEigsSolver.h b/include/Spectra/GenEigsSolver.h index 28b58ea..fcaf4a7 100644 --- a/include/Spectra/GenEigsSolver.h +++ b/include/Spectra/GenEigsSolver.h @@ -120,10 +120,13 @@ class GenEigsSolver : public GenEigsBase { private: using Index = Eigen::Index; + using Scalar = typename OpType::Scalar; + using Complex = std::complex; + using ComplexVector = Eigen::Matrix; public: /// - /// Constructor to create a solver object. + /// 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$: @@ -138,9 +141,10 @@ class GenEigsSolver : public GenEigsBase /// 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) : - GenEigsBase(op, IdentityBOp(), nev, ncv) + GenEigsSolver(OpType& op, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + GenEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) {} }; diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h new file mode 100644 index 0000000..d66924f --- /dev/null +++ b/include/Spectra/LoggerBase.h @@ -0,0 +1,47 @@ + +// Copyright (C) 2018-2021 Yixuan Qiu +// +// 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_LOGGER_BASE_H +#define SPECTRA_LOGGER_BASE_H + +#include +#include // std::vector + +#include "Util/Version.h" +#include "Util/TypeTraits.h" + +namespace Spectra { + +template +class LoggerBase +{ +private: + using Scalar = typename OpType::Scalar; + using Index = Eigen::Index; + using Matrix = Eigen::Matrix; + using Array = Eigen::Array; + using BoolArray = Eigen::Array; + +public: + LoggerBase() {} + + /// + /// Virtual destructor + /// + virtual ~LoggerBase() {} + + // I am not sure what should be in the Data, probably iteration count, residues, current eigenvalues, for davidson maybe subspace size, number_of_converged eigenvalues(). + + /// + /// Virtual logging function + /// + virtual void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged); +}; + +} // namespace Spectra + +#endif // SPECTRA_LOGGER_BASE_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 38db671..778f114 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,6 +11,7 @@ list(APPEND test_target_sources GenEigs.cpp GenEigsRealShift.cpp GenEigsComplexShift.cpp + LoggingGenEigs.cpp Orthogonalization.cpp JDSymEigsBase.cpp JDSymEigsDPRConstructor.cpp diff --git a/test/LoggingGenEigs.cpp b/test/LoggingGenEigs.cpp new file mode 100644 index 0000000..29a2cff --- /dev/null +++ b/test/LoggingGenEigs.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include +#include +#include // Requires C++ 11 +#include // Requires C++ 11 + +#include +#include +#include +#include + +using namespace Spectra; + +#include "catch.hpp" + +using Index = Eigen::Index; +using Matrix = Eigen::MatrixXd; +using Vector = Eigen::VectorXd; +using ComplexMatrix = Eigen::MatrixXcd; +using ComplexVector = Eigen::VectorXcd; +using SpMatrix = Eigen::SparseMatrix; +using BoolArray = Eigen::Array; + +// Derived Logger +template +class DerivedLogger : public LoggerBase +{ + // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. +public: + DerivedLogger(){}; + void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + { + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + std::cout << " Iteration : " << iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; + std::cout << " Size of Krylov subspace : " << subspace_size << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + REQUIRE(residues.size() == current_eigenvalues.size()); + REQUIRE(residues.size() == current_eig_converged.size()); + + std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + for (int i = 0; i < current_eigenvalues.size(); i++) + { + std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + } + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + } +}; + +// Generate random sparse matrix +SpMatrix gen_sparse_data(int n, double prob = 0.5) +{ + SpMatrix mat(n, n); + std::default_random_engine gen; + gen.seed(0); + std::uniform_real_distribution distr(0.0, 1.0); + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + if (distr(gen) < prob) + mat.insert(i, j) = distr(gen) - 0.5; + } + } + return mat; +} + +template +void run_test(const MatType& mat, Solver& eigs, SortRule selection, bool allow_fail = false) +{ + eigs.init(); + // maxit = 300 to reduce running time for failed cases + int nconv = eigs.compute(selection, 300); + int niter = eigs.num_iterations(); + int nops = eigs.num_operations(); + + if (allow_fail && eigs.info() != CompInfo::Successful) + { + WARN("FAILED on this test"); + std::cout << "nconv = " << nconv << std::endl; + std::cout << "niter = " << niter << std::endl; + std::cout << "nops = " << nops << std::endl; + return; + } + else + { + INFO("nconv = " << nconv); + INFO("niter = " << niter); + INFO("nops = " << nops); + REQUIRE(eigs.info() == CompInfo::Successful); + } + + ComplexVector evals = eigs.eigenvalues(); + ComplexMatrix evecs = eigs.eigenvectors(); + + ComplexMatrix resid = mat * evecs - evecs * evals.asDiagonal(); + const double err = resid.array().abs().maxCoeff(); + + INFO("||AU - UD||_inf = " << err); + REQUIRE(err == Approx(0.0).margin(1e-9)); +} + +template +void run_test_sets(const MatType& A, int k, int m) +{ + using OpType = DenseGenMatProd; + + OpType op(A); + std::unique_ptr> logger(new DerivedLogger()); + + GenEigsSolver eigs(op, k, m, std::move(logger)); + run_test(A, eigs, SortRule::LargestMagn); +} + +TEST_CASE("Eigensolver of general real matrix [10x10]", "[eigs_gen]") +{ + std::srand(123); + + const Matrix A = Matrix::Random(10, 10); + int k = 3; + int m = 6; + + run_test_sets(A, k, m); +} From fe4c5d2c0260ccb499899ae83b8f2cb3c2e4fbb2 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Tue, 5 Oct 2021 18:01:40 -0400 Subject: [PATCH 02/16] SymEigs* logging --- include/Spectra/SymEigsBase.h | 28 +++++++++++++++++++++++---- include/Spectra/SymEigsShiftSolver.h | 5 +++-- include/Spectra/SymEigsSolver.h | 6 ++++-- include/Spectra/SymGEigsShiftSolver.h | 15 ++++++++------ include/Spectra/SymGEigsSolver.h | 11 ++++++----- 5 files changed, 46 insertions(+), 19 deletions(-) diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index b6648db..d1bad16 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -12,6 +12,7 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument +#include // std::unique_ptr #include // std::move #include "Util/Version.h" @@ -23,6 +24,7 @@ #include "LinAlg/UpperHessenbergQR.h" #include "LinAlg/TridiagEigen.h" #include "LinAlg/Lanczos.h" +#include "LoggerBase.h" namespace Spectra { @@ -79,8 +81,10 @@ class SymEigsBase private: Matrix m_ritz_vec; // Ritz vectors Vector m_ritz_est; // last row of m_ritz_vec, also called the Ritz estimates + Vector m_resid; BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation + std::unique_ptr> m_logger; // clang-format on // Move rvalue object to the container @@ -150,9 +154,9 @@ class SymEigsBase // 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.array() < thresh); return m_ritz_conv.count(); } @@ -237,7 +241,7 @@ class SymEigsBase /// \cond // If op is an lvalue - SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv) : + SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : m_op(op), m_n(op.rows()), m_nev(nev), @@ -252,10 +256,12 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); + + set_logger(logger); } // If op is an rvalue - SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv) : + SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : m_op_container(create_op_container(std::move(op))), m_op(m_op_container.front()), m_n(m_op.rows()), @@ -271,8 +277,18 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); + + set_logger(logger); } + /// + /// Sets the logger unique_ptr with a user constructed logger object. + /// + void set_logger(std::unique_ptr>& logger) + { + if (logger) + m_logger = std::move(logger); + } /// /// Virtual destructor /// @@ -296,11 +312,13 @@ class SymEigsBase 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; @@ -357,6 +375,8 @@ class SymEigsBase for (i = 0; i < maxit; i++) { nconv = num_converged(tol); + if (m_logger) + m_logger->iteration_log(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); if (nconv >= m_nev) break; diff --git a/include/Spectra/SymEigsShiftSolver.h b/include/Spectra/SymEigsShiftSolver.h index 90fd86c..c483fc7 100644 --- a/include/Spectra/SymEigsShiftSolver.h +++ b/include/Spectra/SymEigsShiftSolver.h @@ -152,6 +152,7 @@ class SymEigsShiftSolver : public SymEigsBase using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Array = Eigen::Array; + using Vector = Eigen::Matrix; using Base = SymEigsBase; using Base::m_nev; @@ -187,8 +188,8 @@ class SymEigsShiftSolver : public SymEigsBase /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) : - Base(op, IdentityBOp(), nev, ncv), + SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigma(sigma) { op.set_shift(m_sigma); diff --git a/include/Spectra/SymEigsSolver.h b/include/Spectra/SymEigsSolver.h index 323c172..5363ed2 100644 --- a/include/Spectra/SymEigsSolver.h +++ b/include/Spectra/SymEigsSolver.h @@ -134,7 +134,9 @@ template > class SymEigsSolver : public SymEigsBase { private: + using Scalar = typename OpType::Scalar; using Index = Eigen::Index; + using Vector = Eigen::Matrix; public: /// @@ -154,8 +156,8 @@ class SymEigsSolver : public SymEigsBase /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymEigsSolver(OpType& op, Index nev, Index ncv) : - SymEigsBase(op, IdentityBOp(), nev, ncv) + SymEigsSolver(OpType& op, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) {} }; diff --git a/include/Spectra/SymGEigsShiftSolver.h b/include/Spectra/SymGEigsShiftSolver.h index 72c9774..8cf03cf 100644 --- a/include/Spectra/SymGEigsShiftSolver.h +++ b/include/Spectra/SymGEigsShiftSolver.h @@ -151,6 +151,7 @@ class SymGEigsShiftSolver : using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Array = Eigen::Array; + using Vector = Eigen::Matrix; using ModeMatOp = SymGEigsShiftInvertOp; using Base = SymEigsBase; @@ -197,8 +198,8 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : - Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} }; @@ -311,6 +312,7 @@ class SymGEigsShiftSolver : using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Array = Eigen::Array; + using Vector = Eigen::Matrix; using ModeMatOp = SymGEigsBucklingOp; using Base = SymEigsBase; @@ -360,8 +362,8 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : - Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} }; @@ -403,6 +405,7 @@ class SymGEigsShiftSolver : using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Array = Eigen::Array; + using Vector = Eigen::Matrix; using ModeMatOp = SymGEigsCayleyOp; using Base = SymEigsBase; @@ -452,8 +455,8 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : - Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} }; diff --git a/include/Spectra/SymGEigsSolver.h b/include/Spectra/SymGEigsSolver.h index ee4cc68..7811345 100644 --- a/include/Spectra/SymGEigsSolver.h +++ b/include/Spectra/SymGEigsSolver.h @@ -185,8 +185,8 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) : - Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv), + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv, std::move(logger)), m_Bop(Bop) {} @@ -253,7 +253,8 @@ class SymGEigsSolver : { private: using Index = Eigen::Index; - + using Scalar = typename OpType::Scalar; + using Vector = Eigen::Matrix; using ModeMatOp = SymGEigsRegInvOp; using Base = SymEigsBase; @@ -280,8 +281,8 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) : - Base(ModeMatOp(op, Bop), Bop, nev, ncv) + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + Base(ModeMatOp(op, Bop), Bop, nev, ncv, std::move(logger)) {} }; From eb940b0c9861b08ba215ad203f0626b3fa5dad09 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Tue, 5 Oct 2021 18:03:17 -0400 Subject: [PATCH 03/16] JDEig* logging --- include/Spectra/DavidsonSymEigsSolver.h | 8 ++++---- include/Spectra/JDSymEigsBase.h | 19 ++++++++++++++++++- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/include/Spectra/DavidsonSymEigsSolver.h b/include/Spectra/DavidsonSymEigsSolver.h index fe9f56c..c97c053 100644 --- a/include/Spectra/DavidsonSymEigsSolver.h +++ b/include/Spectra/DavidsonSymEigsSolver.h @@ -39,8 +39,8 @@ class DavidsonSymEigsSolver : public JDSymEigsBase Vector m_diagonal; public: - DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max) : - JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max) + DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::unique_ptr> logger = nullptr) : + JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max, std::move(logger)) { m_diagonal.resize(this->m_matrix_operator.rows()); for (Index i = 0; i < op.rows(); i++) @@ -49,8 +49,8 @@ class DavidsonSymEigsSolver : public JDSymEigsBase } } - DavidsonSymEigsSolver(OpType& op, Index nev) : - DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev) {} + DavidsonSymEigsSolver(OpType& op, Index nev, std::unique_ptr> logger = nullptr) : + DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, std::move(logger)) {} /// Create initial search space based on the diagonal /// and the spectrum'target (highest or lowest) diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index dca7d08..cccb034 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -12,12 +12,15 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument +#include // std::unique_ptr +#include // std::move #include #include "Util/SelectionRule.h" #include "Util/CompInfo.h" #include "LinAlg/SearchSpace.h" #include "LinAlg/RitzPairs.h" +#include "LoggerBase.h" namespace Spectra { @@ -51,6 +54,7 @@ class JDSymEigsBase private: CompInfo m_info = CompInfo::NotComputed; // status of the computation + std::unique_ptr> m_logger; void check_argument() const { @@ -73,7 +77,7 @@ class JDSymEigsBase } public: - JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max) : + JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::unique_ptr> logger = nullptr) : m_matrix_operator(op), m_number_eigenvalues(nev), m_max_search_space_size(nvec_max < op.rows() ? nvec_max : 10 * m_number_eigenvalues), @@ -82,6 +86,15 @@ class JDSymEigsBase { check_argument(); initialize(); + set_logger(logger); + } + + /// + /// Sets the logger unique_ptr with a user constructed logger object. + /// + void set_logger(std::unique_ptr>& logger) + { + m_logger = std::move(logger); } JDSymEigsBase(OpType& op, Index nev) : @@ -166,6 +179,10 @@ class JDSymEigsBase m_ritz_pairs.sort(selection); bool converged = m_ritz_pairs.check_convergence(tol, m_number_eigenvalues); + if (m_logger) + { + m_logger->iteration_log(niter_, m_ritz_pairs.converged_eigenvalues().count(), m_search_space.size(), eigenvalues(), m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues), m_ritz_pairs.converged_eigenvalues().head(m_number_eigenvalues)); + } if (converged) { m_info = CompInfo::Successful; From 3fe65ac6ab8b7bb5389322a322ddff73254cb5f6 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Tue, 5 Oct 2021 18:04:03 -0400 Subject: [PATCH 04/16] Fix types --- include/Spectra/LoggerBase.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index d66924f..8d0201b 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -16,11 +16,10 @@ namespace Spectra { -template +template class LoggerBase { private: - using Scalar = typename OpType::Scalar; using Index = Eigen::Index; using Matrix = Eigen::Matrix; using Array = Eigen::Array; From d3622e95552c1d955aface588f50f49f875a7ed4 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Tue, 5 Oct 2021 18:04:27 -0400 Subject: [PATCH 05/16] logging tests --- test/CMakeLists.txt | 2 + test/LoggingDavidsonSymEigs.cpp | 144 ++++++++++++++++++++++++++++++++ test/LoggingGenEigs.cpp | 9 +- test/LoggingSymEigs.cpp | 123 +++++++++++++++++++++++++++ 4 files changed, 274 insertions(+), 4 deletions(-) create mode 100644 test/LoggingDavidsonSymEigs.cpp create mode 100644 test/LoggingSymEigs.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 778f114..1df15b3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,7 +11,9 @@ list(APPEND test_target_sources GenEigs.cpp GenEigsRealShift.cpp GenEigsComplexShift.cpp + LoggingDavidsonSymEigs.cpp LoggingGenEigs.cpp + LoggingSymEigs.cpp Orthogonalization.cpp JDSymEigsBase.cpp JDSymEigsDPRConstructor.cpp diff --git a/test/LoggingDavidsonSymEigs.cpp b/test/LoggingDavidsonSymEigs.cpp new file mode 100644 index 0000000..1a7d9d8 --- /dev/null +++ b/test/LoggingDavidsonSymEigs.cpp @@ -0,0 +1,144 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace Spectra; + +#include "catch.hpp" + +using Index = Eigen::Index; + +template +using Matrix = Eigen::Matrix; + +template +using Vector = Eigen::Matrix; + +template +using SpMatrix = Eigen::SparseMatrix; + +using BoolArray = Eigen::Array; + +// Traits to obtain operation type from matrix type +template +struct OpTypeTrait +{ + using Scalar = typename MatType::Scalar; + using OpType = DenseSymMatProd; +}; +template +struct OpTypeTrait> +{ + using OpType = SparseSymMatProd; +}; + +// Derived Logger +template +class DerivedLogger : public LoggerBase +{ + // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. +public: + DerivedLogger(){}; + void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + { + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + std::cout << " Iteration : " << iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; + std::cout << " Size of subspace : " << subspace_size << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + REQUIRE(residues.size() == current_eigenvalues.size()); + REQUIRE(residues.size() == current_eig_converged.size()); + + std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + for (int i = 0; i < current_eigenvalues.size(); i++) + { + std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + } + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + } +}; + +// Generate data for testing +template +Matrix gen_sym_data_dense(int n) +{ + Matrix mat = 0.03 * Matrix::Random(n, n); + Matrix mat1 = mat + mat.transpose(); + for (Eigen::Index i = 0; i < n; i++) + { + mat1(i, i) += i + 1; + } + return mat1; +} + +template +SpMatrix gen_sym_data_sparse(int n) +{ + // Eigen solver only uses the lower triangle of mat, + // so we don't need to make mat symmetric here. + double prob = 0.5; + SpMatrix mat(n, n); + std::default_random_engine gen; + gen.seed(0); + std::uniform_real_distribution distr(0.0, 1.0); + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + if (distr(gen) < prob) + mat.insert(i, j) = 0.1 * (distr(gen) - 0.5); + if (i == j) + { + mat.coeffRef(i, j) = i + 1; + } + } + } + return mat; +} + +template +void run_test(const MatType& mat, int nev, SortRule selection) +{ + using OpType = typename OpTypeTrait::OpType; + using Scalar = typename OpType::Scalar; + OpType op(mat); + std::unique_ptr>> logger(new DerivedLogger>()); + DavidsonSymEigsSolver eigs(op, nev, std::move(logger)); + int nconv = eigs.compute(selection); + + int niter = eigs.num_iterations(); + REQUIRE(nconv == nev); + INFO("nconv = " << nconv); + INFO("niter = " << niter); + REQUIRE(eigs.info() == CompInfo::Successful); + using T = typename OpType::Scalar; + Vector evals = eigs.eigenvalues(); + Matrix evecs = eigs.eigenvectors(); + + Matrix resid = op * evecs - evecs * evals.asDiagonal(); + const T err = resid.array().abs().maxCoeff(); + + INFO("||AU - UD||_inf = " << err); + REQUIRE(err < 100 * Eigen::NumTraits::dummy_precision()); +} + +template +void run_test_set(const MatType& mat, int k) +{ + run_test(mat, k, SortRule::LargestAlge); +} + +TEMPLATE_TEST_CASE("Davidson Solver of sparse symmetric real matrix [1000x1000]", "", double) +{ + std::srand(123); + int k = 10; + const SpMatrix A = gen_sym_data_sparse>(1000); + run_test_set>(A, k); +} diff --git a/test/LoggingGenEigs.cpp b/test/LoggingGenEigs.cpp index 29a2cff..2a6a691 100644 --- a/test/LoggingGenEigs.cpp +++ b/test/LoggingGenEigs.cpp @@ -24,8 +24,8 @@ using SpMatrix = Eigen::SparseMatrix; using BoolArray = Eigen::Array; // Derived Logger -template -class DerivedLogger : public LoggerBase +template +class DerivedLogger : public LoggerBase { // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: @@ -35,7 +35,7 @@ class DerivedLogger : public LoggerBase std::cout << "--------------------------------------------------------------------------------------------" << std::endl; std::cout << " Iteration : " << iteration << std::endl; std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; - std::cout << " Size of Krylov subspace : " << subspace_size << std::endl; + std::cout << " Size of subspace : " << subspace_size << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; REQUIRE(residues.size() == current_eigenvalues.size()); REQUIRE(residues.size() == current_eig_converged.size()); @@ -107,9 +107,10 @@ template void run_test_sets(const MatType& A, int k, int m) { using OpType = DenseGenMatProd; + using Scalar = typename OpType::Scalar; OpType op(A); - std::unique_ptr> logger(new DerivedLogger()); + std::unique_ptr> logger(new DerivedLogger()); GenEigsSolver eigs(op, k, m, std::move(logger)); run_test(A, eigs, SortRule::LargestMagn); diff --git a/test/LoggingSymEigs.cpp b/test/LoggingSymEigs.cpp new file mode 100644 index 0000000..7664362 --- /dev/null +++ b/test/LoggingSymEigs.cpp @@ -0,0 +1,123 @@ +#include +#include +#include +#include +#include +#include // Requires C++ 11 + +#include +#include +#include + +using namespace Spectra; + +#include "catch.hpp" + +using Index = Eigen::Index; +using Matrix = Eigen::MatrixXd; +using Vector = Eigen::VectorXd; +using SpMatrix = Eigen::SparseMatrix; +using BoolArray = Eigen::Array; + +// Derived Logger +template +class DerivedLogger : public LoggerBase +{ + // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. +public: + DerivedLogger(){}; + void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + { + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + std::cout << " Iteration : " << iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; + std::cout << " Size of subspace : " << subspace_size << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + REQUIRE(residues.size() == current_eigenvalues.size()); + REQUIRE(residues.size() == current_eig_converged.size()); + + std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; + std::cout << " ------------------------------------------------------------------------ " << std::endl; + for (int i = 0; i < current_eigenvalues.size(); i++) + { + std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + } + std::cout << "--------------------------------------------------------------------------------------------" << std::endl; + } +}; + +// Generate data for testing +Matrix gen_dense_data(int n) +{ + const Matrix mat = Eigen::MatrixXd::Random(n, n); + return mat + mat.transpose(); +} + +SpMatrix gen_sparse_data(int n, double prob = 0.5) +{ + // Eigen solver only uses the lower triangle of mat, + // so we don't need to make mat symmetric here. + SpMatrix mat(n, n); + std::default_random_engine gen; + gen.seed(0); + std::uniform_real_distribution distr(0.0, 1.0); + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + if (distr(gen) < prob) + mat.insert(i, j) = distr(gen) - 0.5; + } + } + return mat; +} + +template +void run_test(const MatType& mat, Solver& eigs, SortRule selection) +{ + eigs.init(); + int nconv = eigs.compute(selection); + int niter = eigs.num_iterations(); + int nops = eigs.num_operations(); + + INFO("nconv = " << nconv); + INFO("niter = " << niter); + INFO("nops = " << nops); + REQUIRE(eigs.info() == CompInfo::Successful); + + Vector evals = eigs.eigenvalues(); + Matrix evecs = eigs.eigenvectors(); + + Matrix resid = mat.template selfadjointView() * evecs - evecs * evals.asDiagonal(); + const double err = resid.array().abs().maxCoeff(); + + INFO("||AU - UD||_inf = " << err); + REQUIRE(err == Approx(0.0).margin(1e-9)); +} + +template +void run_test_sets(const MatType& mat, int k, int m) +{ + constexpr bool is_dense = std::is_same::value; + using DenseOp = DenseSymMatProd; + using SparseOp = SparseSymMatProd; + using OpType = typename std::conditional::type; + using Scalar = typename OpType::Scalar; + + OpType op(mat); + std::unique_ptr> logger(new DerivedLogger()); + SymEigsSolver eigs(op, k, m, std::move(logger)); + + run_test(mat, eigs, SortRule::LargestMagn); +} + +TEST_CASE("Eigensolver of symmetric real matrix [10x10]", "[eigs_sym]") +{ + std::srand(123); + + const Matrix A = gen_dense_data(10); + int k = 3; + int m = 6; + + run_test_sets(A, k, m); +} From 6f97d87aced32425c7bff8802f66c5fed403d8c6 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Tue, 5 Oct 2021 18:07:18 -0400 Subject: [PATCH 06/16] remove comment --- include/Spectra/LoggerBase.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index 8d0201b..ef96910 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -33,8 +33,6 @@ class LoggerBase /// virtual ~LoggerBase() {} - // I am not sure what should be in the Data, probably iteration count, residues, current eigenvalues, for davidson maybe subspace size, number_of_converged eigenvalues(). - /// /// Virtual logging function /// From d213e64e4209acc23e077c7da2b1bba7df7fa3ab Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Fri, 8 Oct 2021 10:49:24 -0400 Subject: [PATCH 07/16] unique_ptr -> shared_ptr --- include/Spectra/DavidsonSymEigsSolver.h | 4 ++-- include/Spectra/GenEigsBase.h | 10 +++++----- include/Spectra/GenEigsComplexShiftSolver.h | 2 +- include/Spectra/GenEigsRealShiftSolver.h | 2 +- include/Spectra/GenEigsSolver.h | 2 +- include/Spectra/JDSymEigsBase.h | 10 +++++----- include/Spectra/SymEigsBase.h | 12 ++++++------ include/Spectra/SymEigsShiftSolver.h | 2 +- include/Spectra/SymEigsSolver.h | 2 +- include/Spectra/SymGEigsShiftSolver.h | 6 +++--- include/Spectra/SymGEigsSolver.h | 4 ++-- 11 files changed, 28 insertions(+), 28 deletions(-) diff --git a/include/Spectra/DavidsonSymEigsSolver.h b/include/Spectra/DavidsonSymEigsSolver.h index c97c053..2e0f290 100644 --- a/include/Spectra/DavidsonSymEigsSolver.h +++ b/include/Spectra/DavidsonSymEigsSolver.h @@ -39,7 +39,7 @@ class DavidsonSymEigsSolver : public JDSymEigsBase Vector m_diagonal; public: - DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::unique_ptr> logger = nullptr) : + DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger = nullptr) : JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max, std::move(logger)) { m_diagonal.resize(this->m_matrix_operator.rows()); @@ -49,7 +49,7 @@ class DavidsonSymEigsSolver : public JDSymEigsBase } } - DavidsonSymEigsSolver(OpType& op, Index nev, std::unique_ptr> logger = nullptr) : + DavidsonSymEigsSolver(OpType& op, Index nev, std::shared_ptr> logger = nullptr) : DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, std::move(logger)) {} /// Create initial search space based on the diagonal diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 413c8e6..9356f95 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -13,7 +13,7 @@ #include // std::min, std::copy #include // std::complex, std::conj, std::norm, std::abs #include // std::invalid_argument -#include // std::unique_ptr +#include // std::shared_ptr #include "Util/Version.h" #include "Util/TypeTraits.h" @@ -77,7 +77,7 @@ class GenEigsBase private: BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - std::unique_ptr> m_logger; + std::shared_ptr> m_logger; // clang-format on // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part @@ -326,7 +326,7 @@ class GenEigsBase public: /// \cond - GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : m_op(op), m_n(m_op.rows()), m_nev(nev), @@ -346,9 +346,9 @@ class GenEigsBase } /// - /// Sets the logger unique_ptr with a user constructed logger object. + /// Sets the logger shared_ptr with a user constructed logger object. /// - void set_logger(std::unique_ptr>& logger) + void set_logger(std::shared_ptr>& logger) { if (logger) m_logger = std::move(logger); diff --git a/include/Spectra/GenEigsComplexShiftSolver.h b/include/Spectra/GenEigsComplexShiftSolver.h index 1b1e692..dcd8344 100644 --- a/include/Spectra/GenEigsComplexShiftSolver.h +++ b/include/Spectra/GenEigsComplexShiftSolver.h @@ -148,7 +148,7 @@ class GenEigsComplexShiftSolver : public GenEigsBase /// \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, std::unique_ptr> logger = nullptr) : + GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai, std::shared_ptr> logger = nullptr) : Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigmar(sigmar), m_sigmai(sigmai) { diff --git a/include/Spectra/GenEigsRealShiftSolver.h b/include/Spectra/GenEigsRealShiftSolver.h index 52a19af..a1d9c1e 100644 --- a/include/Spectra/GenEigsRealShiftSolver.h +++ b/include/Spectra/GenEigsRealShiftSolver.h @@ -74,7 +74,7 @@ class GenEigsRealShiftSolver : public GenEigsBase /// \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, std::unique_ptr> logger = nullptr) : + GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigma(sigma) { diff --git a/include/Spectra/GenEigsSolver.h b/include/Spectra/GenEigsSolver.h index fcaf4a7..81bc591 100644 --- a/include/Spectra/GenEigsSolver.h +++ b/include/Spectra/GenEigsSolver.h @@ -143,7 +143,7 @@ class GenEigsSolver : public GenEigsBase /// 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, std::unique_ptr> logger = nullptr) : + GenEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : GenEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) {} }; diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index cccb034..cf7dd0d 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -12,7 +12,7 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument -#include // std::unique_ptr +#include // std::shared_ptr #include // std::move #include @@ -54,7 +54,7 @@ class JDSymEigsBase private: CompInfo m_info = CompInfo::NotComputed; // status of the computation - std::unique_ptr> m_logger; + std::shared_ptr> m_logger; void check_argument() const { @@ -77,7 +77,7 @@ class JDSymEigsBase } public: - JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::unique_ptr> logger = nullptr) : + JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger = nullptr) : m_matrix_operator(op), m_number_eigenvalues(nev), m_max_search_space_size(nvec_max < op.rows() ? nvec_max : 10 * m_number_eigenvalues), @@ -90,9 +90,9 @@ class JDSymEigsBase } /// - /// Sets the logger unique_ptr with a user constructed logger object. + /// Sets the logger shared_ptr with a user constructed logger object. /// - void set_logger(std::unique_ptr>& logger) + void set_logger(std::shared_ptr>& logger) { m_logger = std::move(logger); } diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index d1bad16..fddfed0 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -12,7 +12,7 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument -#include // std::unique_ptr +#include // std::shared_ptr #include // std::move #include "Util/Version.h" @@ -84,7 +84,7 @@ class SymEigsBase Vector m_resid; BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - std::unique_ptr> m_logger; + std::shared_ptr> m_logger; // clang-format on // Move rvalue object to the container @@ -241,7 +241,7 @@ class SymEigsBase /// \cond // If op is an lvalue - SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : m_op(op), m_n(op.rows()), m_nev(nev), @@ -261,7 +261,7 @@ class SymEigsBase } // If op is an rvalue - SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : m_op_container(create_op_container(std::move(op))), m_op(m_op_container.front()), m_n(m_op.rows()), @@ -282,9 +282,9 @@ class SymEigsBase } /// - /// Sets the logger unique_ptr with a user constructed logger object. + /// Sets the logger shared_ptr with a user constructed logger object. /// - void set_logger(std::unique_ptr>& logger) + void set_logger(std::shared_ptr>& logger) { if (logger) m_logger = std::move(logger); diff --git a/include/Spectra/SymEigsShiftSolver.h b/include/Spectra/SymEigsShiftSolver.h index c483fc7..bf7c6f6 100644 --- a/include/Spectra/SymEigsShiftSolver.h +++ b/include/Spectra/SymEigsShiftSolver.h @@ -188,7 +188,7 @@ class SymEigsShiftSolver : public SymEigsBase /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : Base(op, IdentityBOp(), nev, ncv, std::move(logger)), m_sigma(sigma) { diff --git a/include/Spectra/SymEigsSolver.h b/include/Spectra/SymEigsSolver.h index 5363ed2..ec03eb6 100644 --- a/include/Spectra/SymEigsSolver.h +++ b/include/Spectra/SymEigsSolver.h @@ -156,7 +156,7 @@ class SymEigsSolver : public SymEigsBase /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymEigsSolver(OpType& op, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : SymEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) {} }; diff --git a/include/Spectra/SymGEigsShiftSolver.h b/include/Spectra/SymGEigsShiftSolver.h index 8cf03cf..43533b2 100644 --- a/include/Spectra/SymGEigsShiftSolver.h +++ b/include/Spectra/SymGEigsShiftSolver.h @@ -198,7 +198,7 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} @@ -362,7 +362,7 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} @@ -455,7 +455,7 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::unique_ptr> logger = nullptr) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} diff --git a/include/Spectra/SymGEigsSolver.h b/include/Spectra/SymGEigsSolver.h index 7811345..93c02c8 100644 --- a/include/Spectra/SymGEigsSolver.h +++ b/include/Spectra/SymGEigsSolver.h @@ -185,7 +185,7 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv, std::move(logger)), m_Bop(Bop) {} @@ -281,7 +281,7 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::unique_ptr> logger = nullptr) : + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : Base(ModeMatOp(op, Bop), Bop, nev, ncv, std::move(logger)) {} }; From efa471d8eaa44f15e909a053315f433d0eb1a904 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Fri, 8 Oct 2021 12:34:55 -0400 Subject: [PATCH 08/16] Separate constructors --- include/Spectra/DavidsonSymEigsSolver.h | 21 ++++- include/Spectra/GenEigsBase.h | 10 ++- include/Spectra/GenEigsComplexShiftSolver.h | 29 ++++++- include/Spectra/GenEigsRealShiftSolver.h | 28 ++++++- include/Spectra/GenEigsSolver.h | 24 +++++- include/Spectra/JDSymEigsBase.h | 15 +++- include/Spectra/SymEigsBase.h | 17 ++-- include/Spectra/SymEigsShiftSolver.h | 29 ++++++- include/Spectra/SymEigsSolver.h | 25 +++++- include/Spectra/SymGEigsShiftSolver.h | 88 +++++++++++++++++++-- include/Spectra/SymGEigsSolver.h | 63 ++++++++++++++- 11 files changed, 314 insertions(+), 35 deletions(-) diff --git a/include/Spectra/DavidsonSymEigsSolver.h b/include/Spectra/DavidsonSymEigsSolver.h index 2e0f290..d2e55df 100644 --- a/include/Spectra/DavidsonSymEigsSolver.h +++ b/include/Spectra/DavidsonSymEigsSolver.h @@ -39,8 +39,8 @@ class DavidsonSymEigsSolver : public JDSymEigsBase Vector m_diagonal; public: - DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger = nullptr) : - JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max, std::move(logger)) + DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max) : + JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max) { m_diagonal.resize(this->m_matrix_operator.rows()); for (Index i = 0; i < op.rows(); i++) @@ -49,8 +49,21 @@ class DavidsonSymEigsSolver : public JDSymEigsBase } } - DavidsonSymEigsSolver(OpType& op, Index nev, std::shared_ptr> logger = nullptr) : - DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, std::move(logger)) {} + DavidsonSymEigsSolver(OpType& op, Index nev) : + DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev) {} + + DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger) : + JDSymEigsBase, 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, std::shared_ptr> logger) : + DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, logger) {} /// Create initial search space based on the diagonal /// and the spectrum'target (highest or lowest) diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 9356f95..7c4e1ec 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -325,8 +325,7 @@ class GenEigsBase public: /// \cond - - GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : + GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv) : m_op(op), m_n(m_op.rows()), m_nev(nev), @@ -341,7 +340,11 @@ class GenEigsBase if (ncv < nev + 2 || ncv > m_n) 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, std::shared_ptr> logger) : + GenEigsBase(op, Bop, nev, ncv) + { set_logger(logger); } @@ -350,8 +353,7 @@ class GenEigsBase /// void set_logger(std::shared_ptr>& logger) { - if (logger) - m_logger = std::move(logger); + m_logger = logger; } /// /// Virtual destructor diff --git a/include/Spectra/GenEigsComplexShiftSolver.h b/include/Spectra/GenEigsComplexShiftSolver.h index dcd8344..9c780f6 100644 --- a/include/Spectra/GenEigsComplexShiftSolver.h +++ b/include/Spectra/GenEigsComplexShiftSolver.h @@ -128,6 +128,31 @@ class GenEigsComplexShiftSolver : public GenEigsBase } public: + /// + /// Constructor to create a eigen solver object using the shift-and-invert mode. + /// + /// \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. + /// + GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai) : + Base(op, IdentityBOp(), nev, ncv), + m_sigmar(sigmar), m_sigmai(sigmai) + { + op.set_shift(m_sigmar, m_sigmai); + } /// /// Constructor to create a eigen solver object using the shift-and-invert mode with logging. /// @@ -148,8 +173,8 @@ class GenEigsComplexShiftSolver : public GenEigsBase /// \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, std::shared_ptr> logger = nullptr) : - Base(op, IdentityBOp(), nev, ncv, std::move(logger)), + GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai, std::shared_ptr> logger) : + Base(op, IdentityBOp(), nev, ncv, logger), m_sigmar(sigmar), m_sigmai(sigmai) { op.set_shift(m_sigmar, m_sigmai); diff --git a/include/Spectra/GenEigsRealShiftSolver.h b/include/Spectra/GenEigsRealShiftSolver.h index a1d9c1e..76980b9 100644 --- a/include/Spectra/GenEigsRealShiftSolver.h +++ b/include/Spectra/GenEigsRealShiftSolver.h @@ -55,6 +55,30 @@ class GenEigsRealShiftSolver : public GenEigsBase } public: + /// + /// Constructor to create a eigen solver object using the shift-and-invert mode. + /// + /// \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. + /// + GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) : + Base(op, IdentityBOp(), nev, ncv), + m_sigma(sigma) + { + op.set_shift(m_sigma); + } /// /// Constructor to create a eigen solver object using the shift-and-invert mode with logging. /// @@ -74,8 +98,8 @@ class GenEigsRealShiftSolver : public GenEigsBase /// \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, std::shared_ptr> logger = nullptr) : - Base(op, IdentityBOp(), nev, ncv, std::move(logger)), + GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + Base(op, IdentityBOp(), nev, ncv, logger), m_sigma(sigma) { op.set_shift(m_sigma); diff --git a/include/Spectra/GenEigsSolver.h b/include/Spectra/GenEigsSolver.h index 81bc591..a212537 100644 --- a/include/Spectra/GenEigsSolver.h +++ b/include/Spectra/GenEigsSolver.h @@ -125,6 +125,26 @@ class GenEigsSolver : public GenEigsBase using ComplexVector = Eigen::Matrix; public: + /// + /// Constructor to create a solver object. + /// + /// \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$. + /// + GenEigsSolver(OpType& op, Index nev, Index ncv) : + GenEigsBase(op, IdentityBOp(), nev, ncv) + {} /// /// Constructor to create a solver object with logging. /// @@ -143,8 +163,8 @@ class GenEigsSolver : public GenEigsBase /// 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, std::shared_ptr> logger = nullptr) : - GenEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) + GenEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger) : + GenEigsBase(op, IdentityBOp(), nev, ncv, logger) {} }; diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index cf7dd0d..b4e2121 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -77,7 +77,18 @@ class JDSymEigsBase } public: - JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger = nullptr) : + JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max) : + 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(); + } + + JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger) : m_matrix_operator(op), m_number_eigenvalues(nev), m_max_search_space_size(nvec_max < op.rows() ? nvec_max : 10 * m_number_eigenvalues), @@ -94,7 +105,7 @@ class JDSymEigsBase /// void set_logger(std::shared_ptr>& logger) { - m_logger = std::move(logger); + m_logger = logger; } JDSymEigsBase(OpType& op, Index nev) : diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index fddfed0..20db2b9 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -241,7 +241,7 @@ class SymEigsBase /// \cond // If op is an lvalue - SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : + SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv) : m_op(op), m_n(op.rows()), m_nev(nev), @@ -256,12 +256,15 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); - + } + SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsBase(op, Bop, nev, ncv) + { set_logger(logger); } // If op is an rvalue - SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : + SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv) : m_op_container(create_op_container(std::move(op))), m_op(m_op_container.front()), m_n(m_op.rows()), @@ -277,7 +280,10 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); - + } + SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsBase(op, Bop, nev, ncv) + { set_logger(logger); } @@ -286,8 +292,7 @@ class SymEigsBase /// void set_logger(std::shared_ptr>& logger) { - if (logger) - m_logger = std::move(logger); + m_logger = logger; } /// /// Virtual destructor diff --git a/include/Spectra/SymEigsShiftSolver.h b/include/Spectra/SymEigsShiftSolver.h index bf7c6f6..f40ffad 100644 --- a/include/Spectra/SymEigsShiftSolver.h +++ b/include/Spectra/SymEigsShiftSolver.h @@ -188,8 +188,33 @@ class SymEigsShiftSolver : public SymEigsBase /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : - Base(op, IdentityBOp(), nev, ncv, std::move(logger)), + SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) : + Base(op, IdentityBOp(), nev, ncv), + m_sigma(sigma) + { + 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 DenseSymShiftSolve, or + /// define their own that implements all the public members + /// as in DenseSymShiftSolve. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param sigma The value of the shift. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + Base(op, IdentityBOp(), nev, ncv, logger), m_sigma(sigma) { op.set_shift(m_sigma); diff --git a/include/Spectra/SymEigsSolver.h b/include/Spectra/SymEigsSolver.h index ec03eb6..074d982 100644 --- a/include/Spectra/SymEigsSolver.h +++ b/include/Spectra/SymEigsSolver.h @@ -156,8 +156,29 @@ class SymEigsSolver : public SymEigsBase /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : - SymEigsBase(op, IdentityBOp(), nev, ncv, std::move(logger)) + SymEigsSolver(OpType& op, Index nev, Index ncv) : + SymEigsBase(op, IdentityBOp(), nev, ncv) + {} + /// + /// Constructor to create a solver object. + /// + /// \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 DenseSymMatProd, or + /// define their own that implements all the public members + /// as in DenseSymMatProd. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsBase(op, IdentityBOp(), nev, ncv, logger) {} }; diff --git a/include/Spectra/SymGEigsShiftSolver.h b/include/Spectra/SymGEigsShiftSolver.h index 43533b2..042dfc8 100644 --- a/include/Spectra/SymGEigsShiftSolver.h +++ b/include/Spectra/SymGEigsShiftSolver.h @@ -198,8 +198,34 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : - Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + m_sigma(sigma) + {} + /// + /// Constructor to create a solver object with logging. + /// + /// \param op The matrix operation object that computes \f$y=(A-\sigma B)^{-1}v\f$ + /// for any vector \f$v\f$. Users could either create the object from the + /// wrapper class SymShiftInvert, or define their own that implements all + /// the public members as in SymShiftInvert. + /// \param Bop The \f$B\f$ matrix operation object that implements the matrix-vector + /// multiplication \f$Bv\f$. Users could either create the object from the + /// wrapper classes such as DenseSymMatProd and SparseSymMatProd, or + /// define their own that implements all the public member functions + /// as in DenseSymMatProd. \f$B\f$ needs to be positive definite. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param sigma The value of the shift. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, logger), m_sigma(sigma) {} }; @@ -362,7 +388,33 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + m_sigma(sigma) + {} + /// + /// Constructor to create a solver object with logging. + /// + /// \param op The matrix operation object that computes \f$y=(K-\sigma K_G)^{-1}v\f$ + /// for any vector \f$v\f$. Users could either create the object from the + /// wrapper class SymShiftInvert, or define their own that implements all + /// the public members as in SymShiftInvert. + /// \param Bop The \f$K\f$ matrix operation object that implements the matrix-vector + /// multiplication \f$Kv\f$. Users could either create the object from the + /// wrapper classes such as DenseSymMatProd and SparseSymMatProd, or + /// define their own that implements all the public member functions + /// as in DenseSymMatProd. \f$K\f$ needs to be positive definite. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param sigma The value of the shift. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} @@ -455,8 +507,34 @@ class SymGEigsShiftSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param sigma The value of the shift. /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger = nullptr) : - Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv), + m_sigma(sigma) + {} + /// + /// Constructor to create a solver object with logging. + /// + /// \param op The matrix operation object that computes \f$y=(A-\sigma B)^{-1}v\f$ + /// for any vector \f$v\f$. Users could either create the object from the + /// wrapper class SymShiftInvert, or define their own that implements all + /// the public members as in SymShiftInvert. + /// \param Bop The \f$B\f$ matrix operation object that implements the matrix-vector + /// multiplication \f$Bv\f$. Users could either create the object from the + /// wrapper classes such as DenseSymMatProd and SparseSymMatProd, or + /// define their own that implements all the public member functions + /// as in DenseSymMatProd. \f$B\f$ needs to be positive definite. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param sigma The value of the shift. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, logger), m_sigma(sigma) {} }; diff --git a/include/Spectra/SymGEigsSolver.h b/include/Spectra/SymGEigsSolver.h index 93c02c8..9a88dcf 100644 --- a/include/Spectra/SymGEigsSolver.h +++ b/include/Spectra/SymGEigsSolver.h @@ -185,8 +185,37 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : - Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv, std::move(logger)), + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) : + Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv), + m_Bop(Bop) + {} + /// + /// Constructor to create a solver object with logging. + /// + /// \param op The \f$A\f$ 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 classes such as DenseSymMatProd, or + /// define their own that implements all the public members + /// as in DenseSymMatProd. + /// \param Bop The \f$B\f$ matrix operation object that represents a Cholesky decomposition of \f$B\f$. + /// It should implement the lower and upper triangular solving operations: + /// calculating \f$L^{-1}v\f$ and \f$(L')^{-1}v\f$ for any vector + /// \f$v\f$, where \f$LL'=B\f$. Users could either + /// create the object from the wrapper classes such as DenseCholesky, or + /// define their own that implements all the public member functions + /// as in DenseCholesky. \f$B\f$ needs to be positive definite. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv, logger), m_Bop(Bop) {} @@ -281,8 +310,34 @@ class SymGEigsSolver : /// in each iteration. This parameter must satisfy \f$nev < ncv \le n\f$, /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger = nullptr) : - Base(ModeMatOp(op, Bop), Bop, nev, ncv, std::move(logger)) + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv) : + Base(ModeMatOp(op, Bop), Bop, nev, ncv) + {} + /// + /// Constructor to create a solver object with logging. + /// + /// \param op The \f$A\f$ 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 classes such as DenseSymMatProd, or + /// define their own that implements all the public members + /// as in DenseSymMatProd. + /// \param Bop The \f$B\f$ matrix operation object that implements the multiplication operation + /// \f$Bv\f$ and the linear equation solving operation \f$B^{-1}v\f$ for any vector \f$v\f$. + /// Users could either create the object from the wrapper class SparseRegularInverse, or + /// define their own that implements all the public member functions + /// as in SparseRegularInverse. \f$B\f$ needs to be positive definite. + /// \param nev Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-1\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 < ncv \le n\f$, + /// and is advised to take \f$ncv \ge 2\cdot nev\f$. + /// \param logger A logging object that inherits from the base class in LoggerBase.h + /// + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + Base(ModeMatOp(op, Bop), Bop, nev, ncv, logger) {} }; From 07e1cbb3ce5dfa32666d396593dcd20c2ae4497d Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Fri, 8 Oct 2021 16:19:45 -0400 Subject: [PATCH 09/16] logging data class --- include/Spectra/GenEigsBase.h | 5 ++++- include/Spectra/JDSymEigsBase.h | 5 ++++- include/Spectra/LoggerBase.h | 24 ++++++++++++++++++++++-- include/Spectra/SymEigsBase.h | 5 ++++- test/LoggingDavidsonSymEigs.cpp | 16 ++++++++-------- test/LoggingGenEigs.cpp | 16 ++++++++-------- test/LoggingSymEigs.cpp | 16 ++++++++-------- 7 files changed, 58 insertions(+), 29 deletions(-) diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 7c4e1ec..9a09042 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -444,7 +444,10 @@ class GenEigsBase { nconv = num_converged(tol); if (m_logger) - m_logger->iteration_log(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + { + const IterationData data(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + m_logger->iteration_log(data); + } if (nconv >= m_nev) break; diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index b4e2121..e3b364f 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -192,7 +192,10 @@ class JDSymEigsBase bool converged = m_ritz_pairs.check_convergence(tol, m_number_eigenvalues); if (m_logger) { - m_logger->iteration_log(niter_, m_ritz_pairs.converged_eigenvalues().count(), m_search_space.size(), eigenvalues(), m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues), m_ritz_pairs.converged_eigenvalues().head(m_number_eigenvalues)); + const Eigen::Array conv_eig = m_ritz_pairs.converged_eigenvalues().head(m_number_eigenvalues); + const Index num_conv = conv_eig.count(); + const IterationData data(niter_, num_conv, m_search_space.size(), eigenvalues(), m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues), conv_eig); + m_logger->iteration_log(data); } if (converged) { diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index ef96910..6a9ba7f 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -17,7 +17,7 @@ namespace Spectra { template -class LoggerBase +class IterationData { private: using Index = Eigen::Index; @@ -25,6 +25,26 @@ class LoggerBase using Array = Eigen::Array; using BoolArray = Eigen::Array; +public: + const Index& iteration; + const Index& number_of_converged; + const Index& subspace_size; + const Vector& current_eigenvalues; + const Vector& residues; + const BoolArray& current_eig_converged; + IterationData(const Index& it, const Index& num_conv, const Index& sub_size, const Vector& cur_eigv, const Vector& res, const BoolArray& cur_eig_conv) : + iteration(it), + number_of_converged(num_conv), + subspace_size(sub_size), + current_eigenvalues(cur_eigv), + residues(res), + current_eig_converged(cur_eig_conv) + {} +}; + +template +class LoggerBase +{ public: LoggerBase() {} @@ -36,7 +56,7 @@ class LoggerBase /// /// Virtual logging function /// - virtual void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged); + virtual void iteration_log(const IterationData& data); }; } // namespace Spectra diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index 20db2b9..5548ed1 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -381,7 +381,10 @@ class SymEigsBase { nconv = num_converged(tol); if (m_logger) - m_logger->iteration_log(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + { + const IterationData data(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + m_logger->iteration_log(data); + } if (nconv >= m_nev) break; diff --git a/test/LoggingDavidsonSymEigs.cpp b/test/LoggingDavidsonSymEigs.cpp index 1a7d9d8..5c28f5d 100644 --- a/test/LoggingDavidsonSymEigs.cpp +++ b/test/LoggingDavidsonSymEigs.cpp @@ -45,21 +45,21 @@ class DerivedLogger : public LoggerBase // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: DerivedLogger(){}; - void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + void iteration_log(const IterationData& data) override { std::cout << "--------------------------------------------------------------------------------------------" << std::endl; - std::cout << " Iteration : " << iteration << std::endl; - std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; - std::cout << " Size of subspace : " << subspace_size << std::endl; + std::cout << " Iteration : " << data.iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; + std::cout << " Size of subspace : " << data.subspace_size << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - REQUIRE(residues.size() == current_eigenvalues.size()); - REQUIRE(residues.size() == current_eig_converged.size()); + REQUIRE(data.residues.size() == data.current_eigenvalues.size()); + REQUIRE(data.residues.size() == data.current_eig_converged.size()); std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - for (int i = 0; i < current_eigenvalues.size(); i++) + for (int i = 0; i < data.current_eigenvalues.size(); i++) { - std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + std::cout << " " << std::setw(20) << data.current_eigenvalues[i] << std::setw(20) << data.current_eig_converged[i] << std::setw(20) << data.residues[i] << std::endl; } std::cout << "--------------------------------------------------------------------------------------------" << std::endl; } diff --git a/test/LoggingGenEigs.cpp b/test/LoggingGenEigs.cpp index 2a6a691..bee10bb 100644 --- a/test/LoggingGenEigs.cpp +++ b/test/LoggingGenEigs.cpp @@ -30,21 +30,21 @@ class DerivedLogger : public LoggerBase // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: DerivedLogger(){}; - void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + void iteration_log(const IterationData& data) override { std::cout << "--------------------------------------------------------------------------------------------" << std::endl; - std::cout << " Iteration : " << iteration << std::endl; - std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; - std::cout << " Size of subspace : " << subspace_size << std::endl; + std::cout << " Iteration : " << data.iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; + std::cout << " Size of subspace : " << data.subspace_size << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - REQUIRE(residues.size() == current_eigenvalues.size()); - REQUIRE(residues.size() == current_eig_converged.size()); + REQUIRE(data.residues.size() == data.current_eigenvalues.size()); + REQUIRE(data.residues.size() == data.current_eig_converged.size()); std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - for (int i = 0; i < current_eigenvalues.size(); i++) + for (int i = 0; i < data.current_eigenvalues.size(); i++) { - std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + std::cout << " " << std::setw(20) << data.current_eigenvalues[i] << std::setw(20) << data.current_eig_converged[i] << std::setw(20) << data.residues[i] << std::endl; } std::cout << "--------------------------------------------------------------------------------------------" << std::endl; } diff --git a/test/LoggingSymEigs.cpp b/test/LoggingSymEigs.cpp index 7664362..6c81095 100644 --- a/test/LoggingSymEigs.cpp +++ b/test/LoggingSymEigs.cpp @@ -26,21 +26,21 @@ class DerivedLogger : public LoggerBase // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: DerivedLogger(){}; - void iteration_log(const Index& iteration, const Index& number_of_converged, const Index& subspace_size, const Vector& current_eigenvalues, const Vector& residues, const BoolArray& current_eig_converged) override + void iteration_log(const IterationData& data) override { std::cout << "--------------------------------------------------------------------------------------------" << std::endl; - std::cout << " Iteration : " << iteration << std::endl; - std::cout << " Number of converged eigenvalues : " << number_of_converged << std::endl; - std::cout << " Size of subspace : " << subspace_size << std::endl; + std::cout << " Iteration : " << data.iteration << std::endl; + std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; + std::cout << " Size of subspace : " << data.subspace_size << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - REQUIRE(residues.size() == current_eigenvalues.size()); - REQUIRE(residues.size() == current_eig_converged.size()); + REQUIRE(data.residues.size() == data.current_eigenvalues.size()); + REQUIRE(data.residues.size() == data.current_eig_converged.size()); std::cout << " " << std::setw(20) << "Current Eigenvalue" << std::setw(20) << "Converged?" << std::setw(20) << "Residue" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; - for (int i = 0; i < current_eigenvalues.size(); i++) + for (int i = 0; i < data.current_eigenvalues.size(); i++) { - std::cout << " " << std::setw(20) << current_eigenvalues[i] << std::setw(20) << current_eig_converged[i] << std::setw(20) << residues[i] << std::endl; + std::cout << " " << std::setw(20) << data.current_eigenvalues[i] << std::setw(20) << data.current_eig_converged[i] << std::setw(20) << data.residues[i] << std::endl; } std::cout << "--------------------------------------------------------------------------------------------" << std::endl; } From b6ecbe2d37d41ae5a12cf30823b30e587ee45e9f Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Sat, 16 Oct 2021 13:09:07 -0400 Subject: [PATCH 10/16] Use raw ptr and prune includes --- include/Spectra/DavidsonSymEigsSolver.h | 4 ++-- include/Spectra/GenEigsBase.h | 9 ++++----- include/Spectra/GenEigsComplexShiftSolver.h | 2 +- include/Spectra/GenEigsRealShiftSolver.h | 2 +- include/Spectra/GenEigsSolver.h | 2 +- include/Spectra/JDSymEigsBase.h | 9 ++++----- include/Spectra/LoggerBase.h | 4 ---- include/Spectra/SymEigsBase.h | 11 +++++------ include/Spectra/SymEigsShiftSolver.h | 2 +- include/Spectra/SymEigsSolver.h | 2 +- include/Spectra/SymGEigsShiftSolver.h | 6 +++--- include/Spectra/SymGEigsSolver.h | 4 ++-- test/LoggingDavidsonSymEigs.cpp | 4 ++-- test/LoggingGenEigs.cpp | 4 ++-- test/LoggingSymEigs.cpp | 4 ++-- 15 files changed, 31 insertions(+), 38 deletions(-) diff --git a/include/Spectra/DavidsonSymEigsSolver.h b/include/Spectra/DavidsonSymEigsSolver.h index d2e55df..9a10f82 100644 --- a/include/Spectra/DavidsonSymEigsSolver.h +++ b/include/Spectra/DavidsonSymEigsSolver.h @@ -52,7 +52,7 @@ class DavidsonSymEigsSolver : public JDSymEigsBase DavidsonSymEigsSolver(OpType& op, Index nev) : DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev) {} - DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger) : + DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max, LoggerBase* logger) : JDSymEigsBase, OpType>(op, nev, nvec_init, nvec_max, logger) { m_diagonal.resize(this->m_matrix_operator.rows()); @@ -62,7 +62,7 @@ class DavidsonSymEigsSolver : public JDSymEigsBase } } - DavidsonSymEigsSolver(OpType& op, Index nev, std::shared_ptr> logger) : + DavidsonSymEigsSolver(OpType& op, Index nev, LoggerBase* logger) : DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev, logger) {} /// Create initial search space based on the diagonal diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 9a09042..eedf0fe 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -13,7 +13,6 @@ #include // std::min, std::copy #include // std::complex, std::conj, std::norm, std::abs #include // std::invalid_argument -#include // std::shared_ptr #include "Util/Version.h" #include "Util/TypeTraits.h" @@ -77,7 +76,7 @@ class GenEigsBase private: BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - std::shared_ptr> m_logger; + LoggerBase *m_logger; // clang-format on // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part @@ -342,16 +341,16 @@ 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, std::shared_ptr> logger) : + GenEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, LoggerBase* logger) : GenEigsBase(op, Bop, nev, ncv) { set_logger(logger); } /// - /// Sets the logger shared_ptr with a user constructed logger object. + /// Sets the logger ptr with a user constructed logger object. /// - void set_logger(std::shared_ptr>& logger) + void set_logger(LoggerBase* logger) { m_logger = logger; } diff --git a/include/Spectra/GenEigsComplexShiftSolver.h b/include/Spectra/GenEigsComplexShiftSolver.h index 9c780f6..18931f6 100644 --- a/include/Spectra/GenEigsComplexShiftSolver.h +++ b/include/Spectra/GenEigsComplexShiftSolver.h @@ -173,7 +173,7 @@ class GenEigsComplexShiftSolver : public GenEigsBase /// \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, std::shared_ptr> logger) : + GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai, LoggerBase* logger) : Base(op, IdentityBOp(), nev, ncv, logger), m_sigmar(sigmar), m_sigmai(sigmai) { diff --git a/include/Spectra/GenEigsRealShiftSolver.h b/include/Spectra/GenEigsRealShiftSolver.h index 76980b9..11952fd 100644 --- a/include/Spectra/GenEigsRealShiftSolver.h +++ b/include/Spectra/GenEigsRealShiftSolver.h @@ -98,7 +98,7 @@ class GenEigsRealShiftSolver : public GenEigsBase /// \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, std::shared_ptr> logger) : + GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, LoggerBase* logger) : Base(op, IdentityBOp(), nev, ncv, logger), m_sigma(sigma) { diff --git a/include/Spectra/GenEigsSolver.h b/include/Spectra/GenEigsSolver.h index a212537..f7fd6b7 100644 --- a/include/Spectra/GenEigsSolver.h +++ b/include/Spectra/GenEigsSolver.h @@ -163,7 +163,7 @@ class GenEigsSolver : public GenEigsBase /// 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, std::shared_ptr> logger) : + GenEigsSolver(OpType& op, Index nev, Index ncv, LoggerBase* logger) : GenEigsBase(op, IdentityBOp(), nev, ncv, logger) {} }; diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index e3b364f..0b6a4b1 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -12,7 +12,6 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument -#include // std::shared_ptr #include // std::move #include @@ -54,7 +53,7 @@ class JDSymEigsBase private: CompInfo m_info = CompInfo::NotComputed; // status of the computation - std::shared_ptr> m_logger; + LoggerBase* m_logger; void check_argument() const { @@ -88,7 +87,7 @@ class JDSymEigsBase initialize(); } - JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, std::shared_ptr> logger) : + JDSymEigsBase(OpType& op, Index nev, Index nvec_init, Index nvec_max, LoggerBase* logger) : m_matrix_operator(op), m_number_eigenvalues(nev), m_max_search_space_size(nvec_max < op.rows() ? nvec_max : 10 * m_number_eigenvalues), @@ -101,9 +100,9 @@ class JDSymEigsBase } /// - /// Sets the logger shared_ptr with a user constructed logger object. + /// Sets the logger ptr with a user constructed logger object. /// - void set_logger(std::shared_ptr>& logger) + void set_logger(LoggerBase* logger) { m_logger = logger; } diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index 6a9ba7f..b1b8a2e 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -9,10 +9,6 @@ #define SPECTRA_LOGGER_BASE_H #include -#include // std::vector - -#include "Util/Version.h" -#include "Util/TypeTraits.h" namespace Spectra { diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index 5548ed1..3b8fadb 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -12,7 +12,6 @@ #include // std::abs, std::pow #include // std::min #include // std::invalid_argument -#include // std::shared_ptr #include // std::move #include "Util/Version.h" @@ -84,7 +83,7 @@ class SymEigsBase Vector m_resid; BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - std::shared_ptr> m_logger; + LoggerBase *m_logger; // clang-format on // Move rvalue object to the container @@ -257,7 +256,7 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); } - SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsBase(OpType& op, const BOpType& Bop, Index nev, Index ncv, LoggerBase* logger) : SymEigsBase(op, Bop, nev, ncv) { set_logger(logger); @@ -281,16 +280,16 @@ class SymEigsBase if (ncv <= nev || ncv > m_n) throw std::invalid_argument("ncv must satisfy nev < ncv <= n, n is the size of matrix"); } - SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsBase(OpType&& op, const BOpType& Bop, Index nev, Index ncv, LoggerBase* logger) : SymEigsBase(op, Bop, nev, ncv) { set_logger(logger); } /// - /// Sets the logger shared_ptr with a user constructed logger object. + /// Sets the logger ptr with a user constructed logger object. /// - void set_logger(std::shared_ptr>& logger) + void set_logger(LoggerBase* logger) { m_logger = logger; } diff --git a/include/Spectra/SymEigsShiftSolver.h b/include/Spectra/SymEigsShiftSolver.h index f40ffad..85487f9 100644 --- a/include/Spectra/SymEigsShiftSolver.h +++ b/include/Spectra/SymEigsShiftSolver.h @@ -213,7 +213,7 @@ class SymEigsShiftSolver : public SymEigsBase /// \param sigma The value of the shift. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + SymEigsShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma, LoggerBase* logger) : Base(op, IdentityBOp(), nev, ncv, logger), m_sigma(sigma) { diff --git a/include/Spectra/SymEigsSolver.h b/include/Spectra/SymEigsSolver.h index 074d982..480222e 100644 --- a/include/Spectra/SymEigsSolver.h +++ b/include/Spectra/SymEigsSolver.h @@ -177,7 +177,7 @@ class SymEigsSolver : public SymEigsBase /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymEigsSolver(OpType& op, Index nev, Index ncv, std::shared_ptr> logger) : + SymEigsSolver(OpType& op, Index nev, Index ncv, LoggerBase* logger) : SymEigsBase(op, IdentityBOp(), nev, ncv, logger) {} }; diff --git a/include/Spectra/SymGEigsShiftSolver.h b/include/Spectra/SymGEigsShiftSolver.h index 042dfc8..f6a4508 100644 --- a/include/Spectra/SymGEigsShiftSolver.h +++ b/include/Spectra/SymGEigsShiftSolver.h @@ -224,7 +224,7 @@ class SymGEigsShiftSolver : /// \param sigma The value of the shift. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, LoggerBase* logger) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, logger), m_sigma(sigma) {} @@ -414,7 +414,7 @@ class SymGEigsShiftSolver : /// \param sigma The value of the shift. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, LoggerBase* logger) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, std::move(logger)), m_sigma(sigma) {} @@ -533,7 +533,7 @@ class SymGEigsShiftSolver : /// \param sigma The value of the shift. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, std::shared_ptr> logger) : + SymGEigsShiftSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, const Scalar& sigma, LoggerBase* logger) : Base(set_shift_and_move(ModeMatOp(op, Bop), sigma), Bop, nev, ncv, logger), m_sigma(sigma) {} diff --git a/include/Spectra/SymGEigsSolver.h b/include/Spectra/SymGEigsSolver.h index 9a88dcf..71d49ed 100644 --- a/include/Spectra/SymGEigsSolver.h +++ b/include/Spectra/SymGEigsSolver.h @@ -214,7 +214,7 @@ class SymGEigsSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, LoggerBase* logger) : Base(ModeMatOp(op, Bop), IdentityBOp(), nev, ncv, logger), m_Bop(Bop) {} @@ -336,7 +336,7 @@ class SymGEigsSolver : /// and is advised to take \f$ncv \ge 2\cdot nev\f$. /// \param logger A logging object that inherits from the base class in LoggerBase.h /// - SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, std::shared_ptr> logger) : + SymGEigsSolver(OpType& op, BOpType& Bop, Index nev, Index ncv, LoggerBase* logger) : Base(ModeMatOp(op, Bop), Bop, nev, ncv, logger) {} }; diff --git a/test/LoggingDavidsonSymEigs.cpp b/test/LoggingDavidsonSymEigs.cpp index 5c28f5d..6aef8aa 100644 --- a/test/LoggingDavidsonSymEigs.cpp +++ b/test/LoggingDavidsonSymEigs.cpp @@ -109,8 +109,8 @@ void run_test(const MatType& mat, int nev, SortRule selection) using OpType = typename OpTypeTrait::OpType; using Scalar = typename OpType::Scalar; OpType op(mat); - std::unique_ptr>> logger(new DerivedLogger>()); - DavidsonSymEigsSolver eigs(op, nev, std::move(logger)); + LoggerBase>* logger(new DerivedLogger>()); + DavidsonSymEigsSolver eigs(op, nev, logger); int nconv = eigs.compute(selection); int niter = eigs.num_iterations(); diff --git a/test/LoggingGenEigs.cpp b/test/LoggingGenEigs.cpp index bee10bb..b17ba88 100644 --- a/test/LoggingGenEigs.cpp +++ b/test/LoggingGenEigs.cpp @@ -110,9 +110,9 @@ void run_test_sets(const MatType& A, int k, int m) using Scalar = typename OpType::Scalar; OpType op(A); - std::unique_ptr> logger(new DerivedLogger()); + LoggerBase* logger(new DerivedLogger()); - GenEigsSolver eigs(op, k, m, std::move(logger)); + GenEigsSolver eigs(op, k, m, logger); run_test(A, eigs, SortRule::LargestMagn); } diff --git a/test/LoggingSymEigs.cpp b/test/LoggingSymEigs.cpp index 6c81095..90b58de 100644 --- a/test/LoggingSymEigs.cpp +++ b/test/LoggingSymEigs.cpp @@ -105,8 +105,8 @@ void run_test_sets(const MatType& mat, int k, int m) using Scalar = typename OpType::Scalar; OpType op(mat); - std::unique_ptr> logger(new DerivedLogger()); - SymEigsSolver eigs(op, k, m, std::move(logger)); + LoggerBase* logger(new DerivedLogger()); + SymEigsSolver eigs(op, k, m, logger); run_test(mat, eigs, SortRule::LargestMagn); } From 90dcfcb0d0ed9d11450f61b64faea9dd04418b92 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Sat, 16 Oct 2021 16:15:31 -0400 Subject: [PATCH 11/16] initialize ptr as nullptr --- include/Spectra/GenEigsBase.h | 2 +- include/Spectra/JDSymEigsBase.h | 2 +- include/Spectra/SymEigsBase.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index eedf0fe..97ef34f 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -76,7 +76,7 @@ class GenEigsBase private: BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - LoggerBase *m_logger; + LoggerBase *m_logger = nullptr; // clang-format on // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index 0b6a4b1..95285c1 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -53,7 +53,7 @@ class JDSymEigsBase private: CompInfo m_info = CompInfo::NotComputed; // status of the computation - LoggerBase* m_logger; + LoggerBase* m_logger = nullptr; void check_argument() const { diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index 3b8fadb..437fcdf 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -83,7 +83,7 @@ class SymEigsBase Vector m_resid; BoolArray m_ritz_conv; // indicator of the convergence of Ritz values CompInfo m_info; // status of the computation - LoggerBase *m_logger; + LoggerBase *m_logger = nullptr; // clang-format on // Move rvalue object to the container From f180bcd09f09c4a2fe37ce6350f6da58a991d4c8 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Sun, 21 Nov 2021 22:56:00 -0500 Subject: [PATCH 12/16] Minor fixes --- include/Spectra/GenEigsBase.h | 3 ++- include/Spectra/JDSymEigsBase.h | 4 +++- include/Spectra/SymEigsBase.h | 3 ++- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 97ef34f..92195c2 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -444,7 +444,8 @@ class GenEigsBase nconv = num_converged(tol); if (m_logger) { - const IterationData data(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + const ComplexVector eigs = m_ritz_val.head(m_nev); + const IterationData data(i, nconv, m_ncv, eigs, m_resid, m_ritz_conv); m_logger->iteration_log(data); } if (nconv >= m_nev) diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index 95285c1..af7e2a6 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -193,7 +193,9 @@ class JDSymEigsBase { const Eigen::Array conv_eig = m_ritz_pairs.converged_eigenvalues().head(m_number_eigenvalues); const Index num_conv = conv_eig.count(); - const IterationData data(niter_, num_conv, m_search_space.size(), eigenvalues(), m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues), conv_eig); + const Vector evals = eigenvalues(); + const Vector res = m_ritz_pairs.residues().colwise().norm().head(m_number_eigenvalues); + const IterationData data(niter_, num_conv, m_search_space.size(), evals, res, conv_eig); m_logger->iteration_log(data); } if (converged) diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index 437fcdf..db9c81c 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -381,7 +381,8 @@ class SymEigsBase nconv = num_converged(tol); if (m_logger) { - const IterationData data(i, nconv, m_ncv, m_ritz_val.head(m_nev), m_resid, m_ritz_conv); + const Vector eigs = m_ritz_val.head(m_nev); + const IterationData data(i, nconv, m_ncv, eigs, m_resid, m_ritz_conv); m_logger->iteration_log(data); } if (nconv >= m_nev) From 8c993e2fa6ab335ea8fde7df7305db77e6976b2a Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Sun, 21 Nov 2021 22:56:17 -0500 Subject: [PATCH 13/16] Logging function pure virtual --- include/Spectra/LoggerBase.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index b1b8a2e..24a51cf 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -52,7 +52,7 @@ class LoggerBase /// /// Virtual logging function /// - virtual void iteration_log(const IterationData& data); + virtual void iteration_log(const IterationData& data) = 0; }; } // namespace Spectra From e3e0869d00adce9f40a665ac75f4810943a91663 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Wed, 18 May 2022 14:36:27 -0400 Subject: [PATCH 14/16] Fix search space size --- include/Spectra/JDSymEigsBase.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index af7e2a6..8e1aeaa 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -195,7 +195,8 @@ class JDSymEigsBase 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 IterationData data(niter_, num_conv, m_search_space.size(), evals, res, conv_eig); + const Index search_space_size = m_search_space.size(); + const IterationData data(niter_, num_conv, search_space_size, evals, res, conv_eig); m_logger->iteration_log(data); } if (converged) From 7bdec45dade1d5fe8a7cea58e37414e4245d5700 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Mon, 29 Aug 2022 10:48:17 -0400 Subject: [PATCH 15/16] Add a logging function that gets called at the start of the iterations. This commit can be rolled back if desired --- include/Spectra/GenEigsBase.h | 19 ++++++++++++++---- include/Spectra/JDSymEigsBase.h | 34 +++++++++++++++++++++++---------- include/Spectra/LoggerBase.h | 3 ++- include/Spectra/SymEigsBase.h | 20 +++++++++++++++---- test/LoggingDavidsonSymEigs.cpp | 25 +++++++++++++++++++++++- test/LoggingGenEigs.cpp | 25 +++++++++++++++++++++++- test/LoggingSymEigs.cpp | 25 +++++++++++++++++++++++- 7 files changed, 129 insertions(+), 22 deletions(-) diff --git a/include/Spectra/GenEigsBase.h b/include/Spectra/GenEigsBase.h index 92195c2..3d7633a 100644 --- a/include/Spectra/GenEigsBase.h +++ b/include/Spectra/GenEigsBase.h @@ -441,18 +441,29 @@ class GenEigsBase Index i, nconv = 0, nev_adj; for (i = 0; i < maxit; i++) { - nconv = num_converged(tol); if (m_logger) { - const ComplexVector eigs = m_ritz_val.head(m_nev); - const IterationData data(i, nconv, m_ncv, eigs, m_resid, m_ritz_conv); - m_logger->iteration_log(data); + m_logger->call_iteration_start(); } + nconv = num_converged(tol); + const ComplexVector eigs = m_ritz_val.head(m_nev); + const IterationData 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); diff --git a/include/Spectra/JDSymEigsBase.h b/include/Spectra/JDSymEigsBase.h index 8e1aeaa..cfd9294 100644 --- a/include/Spectra/JDSymEigsBase.h +++ b/include/Spectra/JDSymEigsBase.h @@ -171,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) @@ -189,30 +194,39 @@ class JDSymEigsBase m_ritz_pairs.sort(selection); bool converged = m_ritz_pairs.check_convergence(tol, m_number_eigenvalues); - if (m_logger) - { - const Eigen::Array 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 data(niter_, num_conv, search_space_size, evals, res, conv_eig); - m_logger->iteration_log(data); - } + const Eigen::Array 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 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(*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().head(m_number_eigenvalues).sum(); } diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index 24a51cf..d11142f 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -49,10 +49,11 @@ class LoggerBase /// virtual ~LoggerBase() {} + virtual void call_iteration_start(); /// /// Virtual logging function /// - virtual void iteration_log(const IterationData& data) = 0; + virtual void call_iteration_end(const IterationData& data) = 0; }; } // namespace Spectra diff --git a/include/Spectra/SymEigsBase.h b/include/Spectra/SymEigsBase.h index db9c81c..b5b38cf 100644 --- a/include/Spectra/SymEigsBase.h +++ b/include/Spectra/SymEigsBase.h @@ -378,18 +378,30 @@ class SymEigsBase Index i, nconv = 0, nev_adj; for (i = 0; i < maxit; i++) { - nconv = num_converged(tol); if (m_logger) { - const Vector eigs = m_ritz_val.head(m_nev); - const IterationData data(i, nconv, m_ncv, eigs, m_resid, m_ritz_conv); - m_logger->iteration_log(data); + m_logger->call_iteration_start(); } + nconv = num_converged(tol); + const Vector eigs = m_ritz_val.head(m_nev); + const IterationData 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); diff --git a/test/LoggingDavidsonSymEigs.cpp b/test/LoggingDavidsonSymEigs.cpp index 6aef8aa..01ac12d 100644 --- a/test/LoggingDavidsonSymEigs.cpp +++ b/test/LoggingDavidsonSymEigs.cpp @@ -44,13 +44,36 @@ class DerivedLogger : public LoggerBase { // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: + std::chrono::time_point start; + std::chrono::time_point end; DerivedLogger(){}; - void iteration_log(const IterationData& data) override + void call_iteration_start() override { + this->start = std::chrono::high_resolution_clock::now(); + } + void call_iteration_end(const IterationData& data) override + { + this->end = std::chrono::high_resolution_clock::now(); + auto duration = this->end - this->start; + auto h = std::chrono::duration_cast(duration); + duration -= h; + auto m = std::chrono::duration_cast(duration); + duration -= m; + auto s = std::chrono::duration_cast(duration); + duration -= s; + auto ms = std::chrono::duration_cast(duration); + duration -= ms; + auto us = std::chrono::duration_cast(duration); + duration -= us; + auto ns = std::chrono::duration_cast(duration); + duration -= ns; + std::stringstream buffer; std::cout << "--------------------------------------------------------------------------------------------" << std::endl; std::cout << " Iteration : " << data.iteration << std::endl; std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; std::cout << " Size of subspace : " << data.subspace_size << std::endl; + std::cout << " Iteration Time " << std::endl; + std::cout << " " << h.count() << "h:" << m.count() << "m:" << s.count() << "s:" << ms.count() << "ms:" << us.count() << "us:" << ns.count() << "ns" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; REQUIRE(data.residues.size() == data.current_eigenvalues.size()); REQUIRE(data.residues.size() == data.current_eig_converged.size()); diff --git a/test/LoggingGenEigs.cpp b/test/LoggingGenEigs.cpp index b17ba88..d758409 100644 --- a/test/LoggingGenEigs.cpp +++ b/test/LoggingGenEigs.cpp @@ -29,13 +29,36 @@ class DerivedLogger : public LoggerBase { // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: + std::chrono::time_point start; + std::chrono::time_point end; DerivedLogger(){}; - void iteration_log(const IterationData& data) override + void call_iteration_start() override { + this->start = std::chrono::high_resolution_clock::now(); + } + void call_iteration_end(const IterationData& data) override + { + this->end = std::chrono::high_resolution_clock::now(); + auto duration = this->end - this->start; + auto h = std::chrono::duration_cast(duration); + duration -= h; + auto m = std::chrono::duration_cast(duration); + duration -= m; + auto s = std::chrono::duration_cast(duration); + duration -= s; + auto ms = std::chrono::duration_cast(duration); + duration -= ms; + auto us = std::chrono::duration_cast(duration); + duration -= us; + auto ns = std::chrono::duration_cast(duration); + duration -= ns; + std::stringstream buffer; std::cout << "--------------------------------------------------------------------------------------------" << std::endl; std::cout << " Iteration : " << data.iteration << std::endl; std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; std::cout << " Size of subspace : " << data.subspace_size << std::endl; + std::cout << " Iteration Time " << std::endl; + std::cout << " " << h.count() << "h:" << m.count() << "m:" << s.count() << "s:" << ms.count() << "ms:" << us.count() << "us:" << ns.count() << "ns" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; REQUIRE(data.residues.size() == data.current_eigenvalues.size()); REQUIRE(data.residues.size() == data.current_eig_converged.size()); diff --git a/test/LoggingSymEigs.cpp b/test/LoggingSymEigs.cpp index 90b58de..983ef1f 100644 --- a/test/LoggingSymEigs.cpp +++ b/test/LoggingSymEigs.cpp @@ -25,13 +25,36 @@ class DerivedLogger : public LoggerBase { // This derived logging class could have some reference to an ostream or call to another class that wraps ostreams etc. public: + std::chrono::time_point start; + std::chrono::time_point end; DerivedLogger(){}; - void iteration_log(const IterationData& data) override + void call_iteration_start() override { + this->start = std::chrono::high_resolution_clock::now(); + } + void call_iteration_end(const IterationData& data) override + { + this->end = std::chrono::high_resolution_clock::now(); + auto duration = this->end - this->start; + auto h = std::chrono::duration_cast(duration); + duration -= h; + auto m = std::chrono::duration_cast(duration); + duration -= m; + auto s = std::chrono::duration_cast(duration); + duration -= s; + auto ms = std::chrono::duration_cast(duration); + duration -= ms; + auto us = std::chrono::duration_cast(duration); + duration -= us; + auto ns = std::chrono::duration_cast(duration); + duration -= ns; + std::stringstream buffer; std::cout << "--------------------------------------------------------------------------------------------" << std::endl; std::cout << " Iteration : " << data.iteration << std::endl; std::cout << " Number of converged eigenvalues : " << data.number_of_converged << std::endl; std::cout << " Size of subspace : " << data.subspace_size << std::endl; + std::cout << " Iteration Time " << std::endl; + std::cout << " " << h.count() << "h:" << m.count() << "m:" << s.count() << "s:" << ms.count() << "ms:" << us.count() << "us:" << ns.count() << "ns" << std::endl; std::cout << " ------------------------------------------------------------------------ " << std::endl; REQUIRE(data.residues.size() == data.current_eigenvalues.size()); REQUIRE(data.residues.size() == data.current_eig_converged.size()); From 88c3b3048b7f0ce1f8caf75752d50795f4cd2030 Mon Sep 17 00:00:00 2001 From: Shiv Upadhyay Date: Mon, 29 Aug 2022 13:05:50 -0400 Subject: [PATCH 16/16] minor fix --- include/Spectra/LoggerBase.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/Spectra/LoggerBase.h b/include/Spectra/LoggerBase.h index d11142f..4a64a15 100644 --- a/include/Spectra/LoggerBase.h +++ b/include/Spectra/LoggerBase.h @@ -49,7 +49,7 @@ class LoggerBase /// virtual ~LoggerBase() {} - virtual void call_iteration_start(); + virtual void call_iteration_start() = 0; /// /// Virtual logging function ///