Large-Scale Eigenvalue Decomposition and SVD with RSpectra

The RSpectra package provides functions eigs() and eigs_sym() to calculate eigenvalues of general and symmetric matrices respectively. If the matrix is known to be symmetric, eigs_sym() is preferred since it guarantees that the eigenvalues are real.

To obtain eigenvalues of a square matrix A, simply call the eigs() or eigs_sym() function, tell it how many eigenvalues are requested (argument k), and which ones are going to be selected (argument which). By default, which = "LM" means to pick the eigenvalues with the largest magnitude (modulus for complex numbers and absolute value for real numbers).

Below we first generate some random matrices and display some of the eigenvalues and eigenvectors:

Performance

diff -Nru r-cran-rspectra-0.13-1/inst/examples/eigs.R r-cran-rspectra-0.15-0/inst/examples/eigs.R --- r-cran-rspectra-0.13-1/inst/examples/eigs.R 2016-02-28 16:04:09.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/examples/eigs.R 2019-05-30 13:06:54.000000000 +0000 @@ -13,8 +13,14 @@ as(x, "dgCMatrix"), as(x, "dgRMatrix")) # Symmetric matrices -sym1 = list(x + t(x)) -sym2 = list(as(x + t(x), "dsyMatrix")) +y = x + t(x) +sym1 = list(y, + as(y, "dgeMatrix"), + as(y, "dgCMatrix"), + as(y, "dgRMatrix")) +sym2 = list(as(y, "dsyMatrix"), + as(y, "dsCMatrix"), + as(as(y, "dsCMatrix"), "dsRMatrix")) ## Test whether the calculated eigenvalues and eigenvectors satisfy ## A * x = lambda * x diff -Nru r-cran-rspectra-0.13-1/inst/examples/svds.R r-cran-rspectra-0.15-0/inst/examples/svds.R --- r-cran-rspectra-0.13-1/inst/examples/svds.R 2016-02-28 16:04:17.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/examples/svds.R 2019-05-30 16:40:30.000000000 +0000 @@ -15,7 +15,10 @@ as(x, "dgRMatrix")) gent = lapply(gen, t) # Symmetric matrices -sym = list(as(crossprod(x), "dsyMatrix")) +y = crossprod(x) +sym = list(as(y, "dsyMatrix"), + as(y, "dsCMatrix"), + as(as(y, "dsCMatrix"), "dsRMatrix")) ## Test whether the calculated (d, u, v) are consistent with svd() ## Return the largest residual @@ -34,7 +37,7 @@ # "True" values gen0 = svd(x) gen0t = svd(t(x)) -sym0 = svd(crossprod(x)) +sym0 = svd(y) ## Capture test result, including error and warning capture = function(expr, env) diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/ComplexShift_sparseMatrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/ComplexShift_sparseMatrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/ComplexShift_sparseMatrix.h 2017-02-05 17:33:33.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/ComplexShift_sparseMatrix.h 2019-05-30 16:33:11.000000000 +0000 @@ -2,13 +2,14 @@ #define COMPLEXSHIFT_SPARSEMATRIX_H #include +#include "ComplexShift.h" template class ComplexShift_sparseMatrix: public ComplexShift { private: typedef Eigen::SparseMatrix SpMat; - typedef Eigen::MappedSparseMatrix MapSpMat; + typedef Eigen::Map MapSpMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/MatProd_sparseMatrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/MatProd_sparseMatrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/MatProd_sparseMatrix.h 2017-02-05 18:49:33.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/MatProd_sparseMatrix.h 2019-05-30 16:33:41.000000000 +0000 @@ -8,7 +8,8 @@ class MatProd_sparseMatrix: public MatProd { private: - typedef Eigen::MappedSparseMatrix MapSpMat; + typedef Eigen::SparseMatrix SpMat; + typedef Eigen::Map MapSpMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/MatProd_sym_sparseMatrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/MatProd_sym_sparseMatrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/MatProd_sym_sparseMatrix.h 2017-02-05 18:51:05.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/MatProd_sym_sparseMatrix.h 2019-05-30 16:31:04.000000000 +0000 @@ -3,12 +3,14 @@ #include #include "MatProd.h" +#include "SparseMatrixMapping.h" template class MatProd_sym_sparseMatrix: public MatProd { private: - typedef Eigen::MappedSparseMatrix MapSpMat; + typedef Eigen::SparseMatrix SpMat; + typedef Eigen::Map MapSpMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; @@ -19,7 +21,7 @@ public: MatProd_sym_sparseMatrix(SEXP mat_, const int nrow_, const char uplo_ = 'L') : - mat(Rcpp::as(mat_)), + mat(map_sparse(mat_)), n(nrow_), uplo(uplo_) {} diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sparseMatrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sparseMatrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sparseMatrix.h 2017-02-05 18:53:26.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sparseMatrix.h 2019-05-30 16:34:04.000000000 +0000 @@ -2,13 +2,14 @@ #define REALSHIFT_SPARSEMATRIX_H #include +#include "RealShift.h" template class RealShift_sparseMatrix: public RealShift { private: typedef Eigen::SparseMatrix SpMat; - typedef Eigen::MappedSparseMatrix MapSpMat; + typedef Eigen::Map MapSpMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; typedef Eigen::SparseLU< Eigen::SparseMatrix > SpLUSolver; diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sym_matrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sym_matrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sym_matrix.h 2017-02-05 19:02:45.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sym_matrix.h 2019-04-06 01:45:06.000000000 +0000 @@ -2,52 +2,79 @@ #define REALSHIFT_SYM_MATRIX_H #include +#include +#include "RealShift.h" class RealShift_sym_matrix: public RealShift { private: typedef Eigen::MatrixXd Matrix; typedef Eigen::VectorXd Vector; + typedef Eigen::VectorXi IntVector; typedef Eigen::Map MapMat; - typedef Eigen::Map MapConstVec; - typedef Eigen::Map MapVec; - typedef Eigen::LDLT LDLTSolver; + typedef Eigen::Map MapConstMat; + typedef Eigen::Map MapConstVec; + typedef Eigen::Map MapVec; - MapMat mat; const int n; const char uplo; - LDLTSolver solver; + Matrix fac; // store the LDLT factorization result + IntVector perm; // store the LDLT factorization result public: RealShift_sym_matrix(SEXP mat_, const int nrow_, const char uplo_ = 'L') : - mat(REAL(mat_), nrow_, nrow_), - n(nrow_), - uplo(uplo_) - {} + n(nrow_), uplo(uplo_), fac(nrow_, nrow_), perm(nrow_) + { + fac.noalias() = MapConstMat(REAL(mat_), nrow_, nrow_); + } int rows() const { return n; } int cols() const { return n; } + // Use LAPACK to do the LDLT factorization void set_shift(double sigma) { - // Backup the diagonal elements - Vector diag = mat.diagonal(); - mat.diagonal().array() -= sigma; - - if(uplo == 'L') - solver.compute(mat.selfadjointView()); - else - solver.compute(mat.selfadjointView()); + fac.diagonal().array() -= sigma; + + int lwork = -1, info; + double blocksize; - mat.diagonal() = diag; + // Inquire the optimal block size + F77_CALL(dsytrf)(&uplo, &n, + fac.data(), &n, + perm.data(), + &blocksize, &lwork, &info); + + if(info != 0) + Rcpp::stop("RealShift_sym_matrix: factorization failed with the given shift"); + + lwork = int(blocksize); + std::vector work; + work.resize(lwork); + + // LDLT factorization + F77_CALL(dsytrf)(&uplo, &n, + fac.data(), &n, + perm.data(), + &work[0], &lwork, &info); + + if(info != 0) + Rcpp::stop("RealShift_sym_matrix: factorization failed with the given shift"); } // y_out = inv(A - sigma * I) * x_in void perform_op(const double* x_in, double* y_out) { - MapConstVec x(x_in, n); - MapVec y(y_out, n); - y.noalias() = solver.solve(x); + std::copy(x_in, x_in + n, y_out); + + const int nrhs = 1; + int info; + F77_CALL(dsytrs)(&uplo, &n, &nrhs, + fac.data(), &n, + perm.data(), y_out, &n, &info); + + if(info != 0) + Rcpp::stop("RealShift_sym_matrix: input vector has illegal values"); } }; diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sym_sparseMatrix.h r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sym_sparseMatrix.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/RealShift_sym_sparseMatrix.h 2017-02-05 18:55:05.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/RealShift_sym_sparseMatrix.h 2019-05-30 16:31:07.000000000 +0000 @@ -2,13 +2,15 @@ #define REALSHIFT_SYM_SPARSEMATRIX_H #include +#include "RealShift.h" +#include "SparseMatrixMapping.h" template class RealShift_sym_sparseMatrix: public RealShift { private: typedef Eigen::SparseMatrix SpMat; - typedef Eigen::MappedSparseMatrix MapSpMat; + typedef Eigen::Map MapSpMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; typedef Eigen::SimplicialLDLT< Eigen::SparseMatrix > SpLDLSolver; @@ -21,7 +23,7 @@ public: RealShift_sym_sparseMatrix(SEXP mat_, const int nrow_, const char uplo_ = 'L') : - mat(Rcpp::as(mat_)), + mat(map_sparse(mat_)), n(nrow_), uplo(uplo_) {} diff -Nru r-cran-rspectra-0.13-1/inst/include/RMatOp/SparseMatrixMapping.h r-cran-rspectra-0.15-0/inst/include/RMatOp/SparseMatrixMapping.h --- r-cran-rspectra-0.13-1/inst/include/RMatOp/SparseMatrixMapping.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/RMatOp/SparseMatrixMapping.h 2019-05-30 16:30:23.000000000 +0000 @@ -0,0 +1,40 @@ +#ifndef MATPROD_SPARSE_MATRIX_MAPPING_H +#define MATPROD_SPARSE_MATRIX_MAPPING_H + +#include + +// Mapping a dgCMatrix/dsCMatrix/dgRMatrix/dsRMatrix R object to Eigen +// Default is ColMajor +template +inline Eigen::Map< Eigen::SparseMatrix > map_sparse(SEXP mat) +{ + Rcpp::S4 obj(mat); + if(!(obj.is("dgCMatrix") || obj.is("dsCMatrix"))) + throw std::invalid_argument("Need S4 class dgCMatrix or dsCMatrix for a mapped sparse matrix"); + + Rcpp::IntegerVector dim(obj.slot("Dim")), i(obj.slot("i")), p(obj.slot("p")); + Rcpp::NumericVector x(obj.slot("x")); + + return Eigen::Map< Eigen::SparseMatrix >( + dim[0], dim[1], p[dim[1]], p.begin(), i.begin(), x.begin() + ); +} + +// Specialization for RowMajor +template <> +inline Eigen::Map< Eigen::SparseMatrix > map_sparse(SEXP mat) +{ + Rcpp::S4 obj(mat); + if(!(obj.is("dgRMatrix") || obj.is("dsRMatrix"))) + throw std::invalid_argument("Need S4 class dgRMatrix or dsRMatrix for a mapped sparse matrix"); + + Rcpp::IntegerVector dim(obj.slot("Dim")), j(obj.slot("j")), p(obj.slot("p")); + Rcpp::NumericVector x(obj.slot("x")); + + return Eigen::Map< Eigen::SparseMatrix >( + dim[0], dim[1], p[dim[1]], p.begin(), j.begin(), x.begin() + ); +} + + +#endif // MATPROD_SPARSE_MATRIX_MAPPING_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsBase.h r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsBase.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsBase.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsBase.h 2019-06-07 23:15:39.000000000 +0000 @@ -0,0 +1,479 @@ +// Copyright (C) 2018-2019 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 GEN_EIGS_BASE_H +#define GEN_EIGS_BASE_H + +#include +#include // std::vector +#include // std::abs, std::pow, std::sqrt +#include // std::min, std::copy +#include // std::complex, std::conj, std::norm, std::abs +#include // std::invalid_argument + +#include "Util/TypeTraits.h" +#include "Util/SelectionRule.h" +#include "Util/CompInfo.h" +#include "Util/SimpleRandom.h" +#include "MatOp/internal/ArnoldiOp.h" +#include "LinAlg/UpperHessenbergQR.h" +#include "LinAlg/DoubleShiftQR.h" +#include "LinAlg/UpperHessenbergEigen.h" +#include "LinAlg/Arnoldi.h" + +namespace Spectra { + + +/// +/// \ingroup EigenSolver +/// +/// This is the base class for general eigen solvers, mainly for internal use. +/// It is kept here to provide the documentation for member functions of concrete eigen solvers +/// such as GenEigsSolver and GenEigsRealShiftSolver. +/// +template < typename Scalar, + int SelectionRule, + typename OpType, + typename BOpType > +class GenEigsBase +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Matrix; + typedef Eigen::Matrix Vector; + typedef Eigen::Array Array; + typedef Eigen::Array BoolArray; + typedef Eigen::Map MapMat; + typedef Eigen::Map MapVec; + typedef Eigen::Map MapConstVec; + + typedef std::complex Complex; + typedef Eigen::Matrix ComplexMatrix; + typedef Eigen::Matrix ComplexVector; + + typedef ArnoldiOp ArnoldiOpType; + typedef Arnoldi ArnoldiFac; + +protected: + OpType* m_op; // object to conduct matrix operation, + // e.g. matrix-vector product + const Index m_n; // dimension of matrix A + const Index m_nev; // number of eigenvalues requested + const Index m_ncv; // dimension of Krylov subspace in the Arnoldi method + Index m_nmatop; // number of matrix operations called + Index m_niter; // number of restarting iterations + + ArnoldiFac m_fac; // Arnoldi factorization + + ComplexVector m_ritz_val; // Ritz values + ComplexMatrix m_ritz_vec; // Ritz vectors + ComplexVector m_ritz_est; // last row of m_ritz_vec + +private: + BoolArray m_ritz_conv; // indicator of the convergence of Ritz values + int m_info; // status of the computation + + const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow + // ~= 1e-307 for the "double" type + const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type + const Scalar m_eps23; // m_eps^(2/3), used to test the convergence + + // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part + // Complex Ritz values have exact conjugate pairs + // So we use exact tests here + static bool is_complex(const Complex& v) { return v.imag() != Scalar(0); } + static bool is_conj(const Complex& v1, const Complex& v2) { return v1 == Eigen::numext::conj(v2); } + + // Implicitly restarted Arnoldi factorization + void restart(Index k) + { + using std::norm; + + if(k >= m_ncv) + return; + + DoubleShiftQR decomp_ds(m_ncv); + UpperHessenbergQR decomp_hb(m_ncv); + Matrix Q = Matrix::Identity(m_ncv, m_ncv); + + for(Index i = k; i < m_ncv; i++) + { + if(is_complex(m_ritz_val[i]) && is_conj(m_ritz_val[i], m_ritz_val[i + 1])) + { + // H - mu * I = Q1 * R1 + // H <- R1 * Q1 + mu * I = Q1' * H * Q1 + // H - conj(mu) * I = Q2 * R2 + // H <- R2 * Q2 + conj(mu) * I = Q2' * H * Q2 + // + // (H - mu * I) * (H - conj(mu) * I) = Q1 * Q2 * R2 * R1 = Q * R + const Scalar s = Scalar(2) * m_ritz_val[i].real(); + const Scalar t = norm(m_ritz_val[i]); + + decomp_ds.compute(m_fac.matrix_H(), s, t); + + // Q -> Q * Qi + decomp_ds.apply_YQ(Q); + // H -> Q'HQ + // Matrix Q = Matrix::Identity(m_ncv, m_ncv); + // decomp_ds.apply_YQ(Q); + // m_fac_H = Q.transpose() * m_fac_H * Q; + m_fac.compress_H(decomp_ds); + + i++; + } else { + // QR decomposition of H - mu * I, mu is real + decomp_hb.compute(m_fac.matrix_H(), m_ritz_val[i].real()); + + // Q -> Q * Qi + decomp_hb.apply_YQ(Q); + // H -> Q'HQ = RQ + mu * I + m_fac.compress_H(decomp_hb); + } + } + + m_fac.compress_V(Q); + m_fac.factorize_from(k, m_ncv, m_nmatop); + + retrieve_ritzpair(); + } + + // Calculates the number of converged Ritz values + Index num_converged(Scalar tol) + { + // thresh = tol * max(m_eps23, abs(theta)), theta for Ritz value + Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_eps23); + Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac.f_norm(); + // Converged "wanted" Ritz values + m_ritz_conv = (resid < thresh); + + return m_ritz_conv.cast().sum(); + } + + // Returns the adjusted nev for restarting + Index nev_adjusted(Index nconv) + { + using std::abs; + + Index nev_new = m_nev; + for(Index i = m_nev; i < m_ncv; i++) + if(abs(m_ritz_est[i]) < m_near_0) nev_new++; + + // Adjust nev_new, according to dnaup2.f line 660~674 in ARPACK + nev_new += std::min(nconv, (m_ncv - nev_new) / 2); + if(nev_new == 1 && m_ncv >= 6) + nev_new = m_ncv / 2; + else if(nev_new == 1 && m_ncv > 3) + nev_new = 2; + + if(nev_new > m_ncv - 2) + nev_new = m_ncv - 2; + + // Increase nev by one if ritz_val[nev - 1] and + // ritz_val[nev] are conjugate pairs + if(is_complex(m_ritz_val[nev_new - 1]) && + is_conj(m_ritz_val[nev_new - 1], m_ritz_val[nev_new])) + { + nev_new++; + } + + return nev_new; + } + + // Retrieves and sorts Ritz values and Ritz vectors + void retrieve_ritzpair() + { + UpperHessenbergEigen decomp(m_fac.matrix_H()); + const ComplexVector& evals = decomp.eigenvalues(); + ComplexMatrix evecs = decomp.eigenvectors(); + + SortEigenvalue sorting(evals.data(), evals.size()); + std::vector ind = sorting.index(); + + // Copy the Ritz values and vectors to m_ritz_val and m_ritz_vec, respectively + for(Index i = 0; i < m_ncv; i++) + { + m_ritz_val[i] = evals[ind[i]]; + m_ritz_est[i] = evecs(m_ncv - 1, ind[i]); + } + for(Index i = 0; i < m_nev; i++) + { + m_ritz_vec.col(i).noalias() = evecs.col(ind[i]); + } + } + +protected: + // Sorts the first nev Ritz pairs in the specified order + // This is used to return the final results + virtual void sort_ritzpair(int sort_rule) + { + // First make sure that we have a valid index vector + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + std::vector ind = sorting.index(); + + switch(sort_rule) + { + case LARGEST_MAGN: + break; + case LARGEST_REAL: + { + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + ind = sorting.index(); + } + break; + case LARGEST_IMAG: + { + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + ind = sorting.index(); + } + break; + case SMALLEST_MAGN: + { + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + ind = sorting.index(); + } + break; + case SMALLEST_REAL: + { + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + ind = sorting.index(); + } + break; + case SMALLEST_IMAG: + { + SortEigenvalue sorting(m_ritz_val.data(), m_nev); + ind = sorting.index(); + } + break; + default: + throw std::invalid_argument("unsupported sorting rule"); + } + + ComplexVector new_ritz_val(m_ncv); + ComplexMatrix new_ritz_vec(m_ncv, m_nev); + BoolArray new_ritz_conv(m_nev); + + for(Index i = 0; i < m_nev; i++) + { + new_ritz_val[i] = m_ritz_val[ind[i]]; + new_ritz_vec.col(i).noalias() = m_ritz_vec.col(ind[i]); + new_ritz_conv[i] = m_ritz_conv[ind[i]]; + } + + m_ritz_val.swap(new_ritz_val); + m_ritz_vec.swap(new_ritz_vec); + m_ritz_conv.swap(new_ritz_conv); + } + +public: + /// \cond + + GenEigsBase(OpType* op, BOpType* Bop, Index nev, Index ncv) : + m_op(op), + m_n(m_op->rows()), + m_nev(nev), + m_ncv(ncv > m_n ? m_n : ncv), + m_nmatop(0), + m_niter(0), + m_fac(ArnoldiOpType(op, Bop), m_ncv), + m_info(NOT_COMPUTED), + m_near_0(TypeTraits::min() * Scalar(10)), + m_eps(Eigen::NumTraits::epsilon()), + m_eps23(Eigen::numext::pow(m_eps, Scalar(2.0) / 3)) + { + if(nev < 1 || nev > m_n - 2) + throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 2, n is the size of matrix"); + + if(ncv < nev + 2 || ncv > m_n) + throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix"); + } + + /// + /// Virtual destructor + /// + virtual ~GenEigsBase() {} + + /// \endcond + + /// + /// Initializes the solver by providing an initial residual vector. + /// + /// \param init_resid Pointer to the initial residual vector. + /// + /// **Spectra** (and also **ARPACK**) uses an iterative algorithm + /// to find eigenvalues. This function allows the user to provide the initial + /// residual vector. + /// + void init(const Scalar* init_resid) + { + // Reset all matrices/vectors to zero + m_ritz_val.resize(m_ncv); + m_ritz_vec.resize(m_ncv, m_nev); + m_ritz_est.resize(m_ncv); + m_ritz_conv.resize(m_nev); + + m_ritz_val.setZero(); + m_ritz_vec.setZero(); + m_ritz_est.setZero(); + m_ritz_conv.setZero(); + + m_nmatop = 0; + m_niter = 0; + + // Initialize the Arnoldi factorization + MapConstVec v0(init_resid, m_n); + m_fac.init(v0, m_nmatop); + } + + /// + /// Initializes the solver by providing a random initial residual vector. + /// + /// This overloaded function generates a random initial residual vector + /// (with a fixed random seed) for the algorithm. Elements in the vector + /// follow independent Uniform(-0.5, 0.5) distribution. + /// + void init() + { + SimpleRandom rng(0); + Vector init_resid = rng.random_vec(m_n); + init(init_resid.data()); + } + + /// + /// Conducts the major computation procedure. + /// + /// \param maxit Maximum number of iterations allowed in the algorithm. + /// \param tol Precision parameter for the calculated eigenvalues. + /// \param sort_rule Rule to sort the eigenvalues and eigenvectors. + /// Supported values are + /// `Spectra::LARGEST_MAGN`, `Spectra::LARGEST_REAL`, + /// `Spectra::LARGEST_IMAG`, `Spectra::SMALLEST_MAGN`, + /// `Spectra::SMALLEST_REAL` and `Spectra::SMALLEST_IMAG`, + /// for example `LARGEST_MAGN` indicates that eigenvalues + /// with largest magnitude come first. + /// Note that this argument is only used to + /// **sort** the final result, and the **selection** rule + /// (e.g. selecting the largest or smallest eigenvalues in the + /// full spectrum) is specified by the template parameter + /// `SelectionRule` of GenEigsSolver. + /// + /// \return Number of converged eigenvalues. + /// + Index compute(Index maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_MAGN) + { + // The m-step Arnoldi factorization + m_fac.factorize_from(1, m_ncv, m_nmatop); + retrieve_ritzpair(); + // Restarting + Index i, nconv = 0, nev_adj; + for(i = 0; i < maxit; i++) + { + nconv = num_converged(tol); + if(nconv >= m_nev) + break; + + nev_adj = nev_adjusted(nconv); + restart(nev_adj); + } + // Sorting results + sort_ritzpair(sort_rule); + + m_niter += i + 1; + m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING; + + return std::min(m_nev, nconv); + } + + /// + /// Returns the status of the computation. + /// The full list of enumeration values can be found in \ref Enumerations. + /// + int info() const { return m_info; } + + /// + /// Returns the number of iterations used in the computation. + /// + Index num_iterations() const { return m_niter; } + + /// + /// Returns the number of matrix operations used in the computation. + /// + Index num_operations() const { return m_nmatop; } + + /// + /// Returns the converged eigenvalues. + /// + /// \return A complex-valued vector containing the eigenvalues. + /// Returned vector type will be `Eigen::Vector, ...>`, depending on + /// the template parameter `Scalar` defined. + /// + ComplexVector eigenvalues() const + { + const Index nconv = m_ritz_conv.cast().sum(); + ComplexVector res(nconv); + + if(!nconv) + return res; + + Index j = 0; + for(Index i = 0; i < m_nev; i++) + { + if(m_ritz_conv[i]) + { + res[j] = m_ritz_val[i]; + j++; + } + } + + return res; + } + + /// + /// Returns the eigenvectors associated with the converged eigenvalues. + /// + /// \param nvec The number of eigenvectors to return. + /// + /// \return A complex-valued matrix containing the eigenvectors. + /// Returned matrix type will be `Eigen::Matrix, ...>`, + /// depending on the template parameter `Scalar` defined. + /// + ComplexMatrix eigenvectors(Index nvec) const + { + const Index nconv = m_ritz_conv.cast().sum(); + nvec = std::min(nvec, nconv); + ComplexMatrix res(m_n, nvec); + + if(!nvec) + return res; + + ComplexMatrix ritz_vec_conv(m_ncv, nvec); + Index j = 0; + for(Index i = 0; i < m_nev && j < nvec; i++) + { + if(m_ritz_conv[i]) + { + ritz_vec_conv.col(j).noalias() = m_ritz_vec.col(i); + j++; + } + } + + res.noalias() = m_fac.matrix_V() * ritz_vec_conv; + + return res; + } + + /// + /// Returns all converged eigenvectors. + /// + ComplexMatrix eigenvectors() const + { + return eigenvectors(m_nev); + } +}; + + +} // namespace Spectra + +#endif // GEN_EIGS_BASE_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsComplexShiftSolver.h r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsComplexShiftSolver.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsComplexShiftSolver.h 2018-05-22 13:09:59.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsComplexShiftSolver.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -7,7 +7,10 @@ #ifndef GEN_EIGS_COMPLEX_SHIFT_SOLVER_H #define GEN_EIGS_COMPLEX_SHIFT_SOLVER_H -#include "GenEigsSolver.h" +#include + +#include "GenEigsBase.h" +#include "Util/SelectionRule.h" #include "MatOp/DenseGenComplexShiftSolve.h" namespace Spectra { @@ -29,15 +32,16 @@ /// \ref Enumerations. /// \tparam OpType The name of the matrix operation class. Users could either /// use the DenseGenComplexShiftSolve wrapper class, or define their -/// own that impelemnts all the public member functions as in +/// own that implements all the public member functions as in /// DenseGenComplexShiftSolve. /// template > -class GenEigsComplexShiftSolver: public GenEigsSolver +class GenEigsComplexShiftSolver: public GenEigsBase { private: + typedef Eigen::Index Index; typedef std::complex Complex; typedef Eigen::Matrix Vector; typedef Eigen::Matrix ComplexVector; @@ -45,14 +49,14 @@ const Scalar m_sigmar; const Scalar m_sigmai; - // First transform back the ritz values, and then sort + // First transform back the Ritz values, and then sort void sort_ritzpair(int sort_rule) { using std::abs; using std::sqrt; using std::norm; - // The eigenvalus we get from the iteration is + // The eigenvalues we get from the iteration is // nu = 0.5 * (1 / (lambda - sigma) + 1 / (lambda - conj(sigma))) // So the eigenvalues of the original problem is // 1 \pm sqrt(1 - 4 * nu^2 * sigmai^2) @@ -79,10 +83,10 @@ // Calculate inv(A - r * I) * vj Vector v_real(this->m_n), v_imag(this->m_n), OPv_real(this->m_n), OPv_imag(this->m_n); const Scalar eps = Eigen::NumTraits::epsilon(); - for(int i = 0; i < this->m_nev; i++) + for(Index i = 0; i < this->m_nev; i++) { - v_real.noalias() = this->m_fac_V * this->m_ritz_vec.col(i).real(); - v_imag.noalias() = this->m_fac_V * this->m_ritz_vec.col(i).imag(); + v_real.noalias() = this->m_fac.matrix_V() * this->m_ritz_vec.col(i).real(); + v_imag.noalias() = this->m_fac.matrix_V() * this->m_ritz_vec.col(i).imag(); this->m_op->perform_op(v_real.data(), OPv_real.data()); this->m_op->perform_op(v_imag.data(), OPv_imag.data()); @@ -116,7 +120,7 @@ } } - GenEigsSolver::sort_ritzpair(sort_rule); + GenEigsBase::sort_ritzpair(sort_rule); } public: /// @@ -126,7 +130,7 @@ /// 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 DenseGenComplexShiftSolve wrapper class, or - /// define their own that impelemnts all the public member functions + /// define their own that implements all the public member functions /// 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. @@ -138,8 +142,8 @@ /// \param sigmar The real part of the shift. /// \param sigmai The imaginary part of the shift. /// - GenEigsComplexShiftSolver(OpType* op, int nev, int ncv, const Scalar& sigmar, const Scalar& sigmai) : - GenEigsSolver(op, nev, ncv), + GenEigsComplexShiftSolver(OpType* op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai) : + GenEigsBase(op, NULL, nev, ncv), m_sigmar(sigmar), m_sigmai(sigmai) { this->m_op->set_shift(m_sigmar, m_sigmai); diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsRealShiftSolver.h r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsRealShiftSolver.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsRealShiftSolver.h 2018-05-22 13:09:59.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsRealShiftSolver.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -7,7 +7,10 @@ #ifndef GEN_EIGS_REAL_SHIFT_SOLVER_H #define GEN_EIGS_REAL_SHIFT_SOLVER_H -#include "GenEigsSolver.h" +#include + +#include "GenEigsBase.h" +#include "Util/SelectionRule.h" #include "MatOp/DenseGenRealShiftSolve.h" namespace Spectra { @@ -30,28 +33,29 @@ /// \tparam OpType The name of the matrix operation class. Users could either /// use the wrapper classes such as DenseGenRealShiftSolve and /// SparseGenRealShiftSolve, or define their -/// own that impelemnts all the public member functions as in +/// own that implements all the public member functions as in /// DenseGenRealShiftSolve. /// template > -class GenEigsRealShiftSolver: public GenEigsSolver +class GenEigsRealShiftSolver: public GenEigsBase { private: + typedef Eigen::Index Index; typedef std::complex Complex; typedef Eigen::Array ComplexArray; const Scalar m_sigma; - // First transform back the ritz values, and then sort + // First transform back the Ritz values, and then sort void sort_ritzpair(int sort_rule) { - // The eigenvalus we get from the iteration is nu = 1 / (lambda - sigma) + // The eigenvalues we get from the iteration is nu = 1 / (lambda - sigma) // So the eigenvalues of the original problem is lambda = 1 / nu + sigma ComplexArray ritz_val_org = Scalar(1.0) / this->m_ritz_val.head(this->m_nev).array() + m_sigma; this->m_ritz_val.head(this->m_nev) = ritz_val_org; - GenEigsSolver::sort_ritzpair(sort_rule); + GenEigsBase::sort_ritzpair(sort_rule); } public: /// @@ -61,7 +65,7 @@ /// 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 impelemnts all the public member functions + /// define their own that implements all the public member functions /// 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. @@ -72,8 +76,8 @@ /// and is advised to take \f$ncv \ge 2\cdot nev + 1\f$. /// \param sigma The real-valued shift. /// - GenEigsRealShiftSolver(OpType* op, int nev, int ncv, Scalar sigma) : - GenEigsSolver(op, nev, ncv), + GenEigsRealShiftSolver(OpType* op, Index nev, Index ncv, Scalar sigma) : + GenEigsBase(op, NULL, nev, ncv), m_sigma(sigma) { this->m_op->set_shift(m_sigma); diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsSolver.h r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsSolver.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/GenEigsSolver.h 2018-05-22 20:39:13.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/GenEigsSolver.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -8,19 +8,9 @@ #define GEN_EIGS_SOLVER_H #include -#include // std::vector -#include // std::abs, std::pow, std::sqrt -#include // std::min, std::copy -#include // std::complex, std::conj, std::norm, std::abs -#include // std::invalid_argument -#include "Util/TypeTraits.h" +#include "GenEigsBase.h" #include "Util/SelectionRule.h" -#include "Util/CompInfo.h" -#include "Util/SimpleRandom.h" -#include "LinAlg/UpperHessenbergQR.h" -#include "LinAlg/UpperHessenbergEigen.h" -#include "LinAlg/DoubleShiftQR.h" #include "MatOp/DenseGenMatProd.h" namespace Spectra { @@ -46,14 +36,15 @@ /// \tparam OpType The name of the matrix operation class. Users could either /// use the wrapper classes such as DenseGenMatProd and /// SparseGenMatProd, or define their -/// own that impelemnts all the public member functions as in +/// own that implements all the public member functions as in /// DenseGenMatProd. /// /// An example that illustrates the usage of GenEigsSolver is give below: /// /// \code{.cpp} /// #include -/// #include // Also includes +/// #include +/// // is implicitly included /// #include /// /// using namespace Spectra; @@ -90,8 +81,8 @@ /// \code{.cpp} /// #include /// #include -/// #include -/// #include +/// #include +/// #include /// #include /// /// using namespace Spectra; @@ -135,358 +126,10 @@ template < typename Scalar = double, int SelectionRule = LARGEST_MAGN, typename OpType = DenseGenMatProd > -class GenEigsSolver +class GenEigsSolver: public GenEigsBase { private: - typedef Eigen::Matrix Matrix; - typedef Eigen::Matrix Vector; - typedef Eigen::Array Array; - typedef Eigen::Array BoolArray; - typedef Eigen::Map MapMat; - typedef Eigen::Map MapVec; - - typedef std::complex Complex; - typedef Eigen::Matrix ComplexMatrix; - typedef Eigen::Matrix ComplexVector; - -protected: - OpType* m_op; // object to conduct matrix operation, - // e.g. matrix-vector product - const int m_n; // dimension of matrix A - const int m_nev; // number of eigenvalues requested - const int m_ncv; // dimension of Krylov subspace in the Arnoldi method - int m_nmatop; // number of matrix operations called - int m_niter; // number of restarting iterations - - Matrix m_fac_V; // V matrix in the Arnoldi factorization - Matrix m_fac_H; // H matrix in the Arnoldi factorization - Vector m_fac_f; // residual in the Arnoldi factorization - - ComplexVector m_ritz_val; // Ritz values - ComplexMatrix m_ritz_vec; // Ritz vectors - ComplexVector m_ritz_est; // last row of m_ritz_vec - -private: - BoolArray m_ritz_conv; // indicator of the convergence of Ritz values - int m_info; // status of the computation - - const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow - // ~= 1e-307 for the "double" type - const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type - const Scalar m_eps23; // m_eps^(2/3), used to test the convergence - - // Given orthonormal basis functions V, find a nonzero vector f such that V'f = 0 - // Assume that f has been properly allocated - void expand_basis(const MapMat& V, const int seed, Vector& f, Scalar& fnorm) - { - using std::sqrt; - - const Scalar thresh = m_eps * sqrt(Scalar(m_n)); - for(int iter = 0; iter < 5; iter++) - { - // Randomly generate a new vector and orthogonalize it against V - SimpleRandom rng(seed + 123 * iter); - f.noalias() = rng.random_vec(m_n); - // f <- f - V * V' * f, so that f is orthogonal to V - Vector Vf = V.transpose() * f; - f -= V * Vf; - // fnorm <- ||f|| - fnorm = m_fac_f.norm(); - - // If fnorm is too close to zero, we try a new random vector, - // otherwise return the result - if(fnorm >= thresh) - return; - } - } - - // Arnoldi factorization starting from step-k - void factorize_from(int from_k, int to_m, const Vector& fk) - { - using std::sqrt; - - if(to_m <= from_k) return; - - const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n)); - m_fac_f.noalias() = fk; - - // Pre-allocate Vf - Vector Vf(to_m); - Vector w(m_n); - Scalar beta = m_fac_f.norm(); - // Keep the upperleft k x k submatrix of H and set other elements to 0 - m_fac_H.rightCols(m_ncv - from_k).setZero(); - m_fac_H.block(from_k, 0, m_ncv - from_k, from_k).setZero(); - for(int i = from_k; i <= to_m - 1; i++) - { - bool restart = false; - // If beta = 0, then the next V is not full rank - // We need to generate a new residual vector that is orthogonal - // to the current V, which we call a restart - if(beta < m_near_0) - { - MapMat V(m_fac_V.data(), m_n, i); // The first i columns - expand_basis(V, 2 * i, m_fac_f, beta); - restart = true; - } - - // v <- f / ||f|| - m_fac_V.col(i).noalias() = m_fac_f / beta; // The (i+1)-th column - - // Note that H[i+1, i] equals to the unrestarted beta - m_fac_H(i, i - 1) = restart ? Scalar(0) : beta; - - // w <- A * v, v = m_fac_V.col(i) - m_op->perform_op(&m_fac_V(0, i), w.data()); - m_nmatop++; - - const int i1 = i + 1; - // First i+1 columns of V - MapMat Vs(m_fac_V.data(), m_n, i1); - // h = m_fac_H(0:i, i) - MapVec h(&m_fac_H(0, i), i1); - // h <- V' * w - h.noalias() = Vs.transpose() * w; - - // f <- w - V * h - m_fac_f.noalias() = w - Vs * h; - beta = m_fac_f.norm(); - - if(beta > Scalar(0.717) * h.norm()) - continue; - - // f/||f|| is going to be the next column of V, so we need to test - // whether V' * (f/||f||) ~= 0 - Vf.head(i1).noalias() = Vs.transpose() * m_fac_f; - Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); - // If not, iteratively correct the residual - int count = 0; - while(count < 5 && ortho_err > m_eps * beta) - { - // There is an edge case: when beta=||f|| is close to zero, f mostly consists - // of noises of rounding errors, so the test [ortho_err < eps * beta] is very - // likely to fail. In particular, if beta=0, then the test is ensured to fail. - // Hence when this happens, we force f to be zero, and then restart in the - // next iteration. - if(beta < beta_thresh) - { - m_fac_f.setZero(); - beta = Scalar(0); - break; - } - - // f <- f - V * Vf - m_fac_f.noalias() -= Vs * Vf.head(i1); - // h <- h + Vf - h.noalias() += Vf.head(i1); - // beta <- ||f|| - beta = m_fac_f.norm(); - - Vf.head(i1).noalias() = Vs.transpose() * m_fac_f; - ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); - count++; - } - } - } - - // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part - // Complex Ritz values have exact conjugate pairs - // So we use exact tests here - static bool is_complex(const Complex& v) { return v.imag() != Scalar(0); } - static bool is_conj(const Complex& v1, const Complex& v2) { return v1 == Eigen::numext::conj(v2); } - - // Implicitly restarted Arnoldi factorization - void restart(int k) - { - using std::norm; - - if(k >= m_ncv) - return; - - DoubleShiftQR decomp_ds(m_ncv); - UpperHessenbergQR decomp_hb; - Matrix Q = Matrix::Identity(m_ncv, m_ncv); - - for(int i = k; i < m_ncv; i++) - { - if(is_complex(m_ritz_val[i]) && is_conj(m_ritz_val[i], m_ritz_val[i + 1])) - { - // H - mu * I = Q1 * R1 - // H <- R1 * Q1 + mu * I = Q1' * H * Q1 - // H - conj(mu) * I = Q2 * R2 - // H <- R2 * Q2 + conj(mu) * I = Q2' * H * Q2 - // - // (H - mu * I) * (H - conj(mu) * I) = Q1 * Q2 * R2 * R1 = Q * R - const Scalar s = Scalar(2) * m_ritz_val[i].real(); - const Scalar t = norm(m_ritz_val[i]); - - decomp_ds.compute(m_fac_H, s, t); - - // Q -> Q * Qi - decomp_ds.apply_YQ(Q); - // H -> Q'HQ - // Matrix Q = Matrix::Identity(m_ncv, m_ncv); - // decomp_ds.apply_YQ(Q); - // m_fac_H = Q.transpose() * m_fac_H * Q; - m_fac_H.noalias() = decomp_ds.matrix_QtHQ(); - - i++; - } else { - // QR decomposition of H - mu * I, mu is real - m_fac_H.diagonal().array() -= m_ritz_val[i].real(); - decomp_hb.compute(m_fac_H); - - // Q -> Q * Qi - decomp_hb.apply_YQ(Q); - // H -> Q'HQ = RQ + mu * I - decomp_hb.matrix_RQ(m_fac_H); - m_fac_H.diagonal().array() += m_ritz_val[i].real(); - } - } - // V -> VQ, only need to update the first k+1 columns - // Q has some elements being zero - // The first (ncv - k + i) elements of the i-th column of Q are non-zero - Matrix Vs(m_n, k + 1); - for(int i = 0; i < k; i++) - { - const int nnz = m_ncv - k + i + 1; - MapVec q(&Q(0, i), nnz); - Vs.col(i).noalias() = m_fac_V.leftCols(nnz) * q; - } - Vs.col(k).noalias() = m_fac_V * Q.col(k); - m_fac_V.leftCols(k + 1).noalias() = Vs; - - const Vector fk = m_fac_f * Q(m_ncv - 1, k - 1) + m_fac_V.col(k) * m_fac_H(k, k - 1); - factorize_from(k, m_ncv, fk); - retrieve_ritzpair(); - } - - // Calculates the number of converged Ritz values - int num_converged(Scalar tol) - { - // thresh = tol * max(m_eps23, abs(theta)), theta for Ritz value - Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_eps23); - Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac_f.norm(); - // Converged "wanted" Ritz values - m_ritz_conv = (resid < thresh); - - return m_ritz_conv.cast().sum(); - } - - // Returns the adjusted nev for restarting - int nev_adjusted(int nconv) - { - using std::abs; - - int nev_new = m_nev; - for(int i = m_nev; i < m_ncv; i++) - if(abs(m_ritz_est[i]) < m_near_0) nev_new++; - - // Adjust nev_new, according to dnaup2.f line 660~674 in ARPACK - nev_new += std::min(nconv, (m_ncv - nev_new) / 2); - if(nev_new == 1 && m_ncv >= 6) - nev_new = m_ncv / 2; - else if(nev_new == 1 && m_ncv > 3) - nev_new = 2; - - if(nev_new > m_ncv - 2) - nev_new = m_ncv - 2; - - // Increase nev by one if ritz_val[nev - 1] and - // ritz_val[nev] are conjugate pairs - if(is_complex(m_ritz_val[nev_new - 1]) && - is_conj(m_ritz_val[nev_new - 1], m_ritz_val[nev_new])) - { - nev_new++; - } - - return nev_new; - } - - // Retrieves and sorts Ritz values and Ritz vectors - void retrieve_ritzpair() - { - UpperHessenbergEigen decomp(m_fac_H); - const ComplexVector& evals = decomp.eigenvalues(); - ComplexMatrix evecs = decomp.eigenvectors(); - - SortEigenvalue sorting(evals.data(), evals.size()); - std::vector ind = sorting.index(); - - // Copy the Ritz values and vectors to m_ritz_val and m_ritz_vec, respectively - for(int i = 0; i < m_ncv; i++) - { - m_ritz_val[i] = evals[ind[i]]; - m_ritz_est[i] = evecs(m_ncv - 1, ind[i]); - } - for(int i = 0; i < m_nev; i++) - { - m_ritz_vec.col(i).noalias() = evecs.col(ind[i]); - } - } - -protected: - // Sorts the first nev Ritz pairs in the specified order - // This is used to return the final results - virtual void sort_ritzpair(int sort_rule) - { - // First make sure that we have a valid index vector - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - std::vector ind = sorting.index(); - - switch(sort_rule) - { - case LARGEST_MAGN: - break; - case LARGEST_REAL: - { - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - ind = sorting.index(); - } - break; - case LARGEST_IMAG: - { - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - ind = sorting.index(); - } - break; - case SMALLEST_MAGN: - { - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - ind = sorting.index(); - } - break; - case SMALLEST_REAL: - { - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - ind = sorting.index(); - } - break; - case SMALLEST_IMAG: - { - SortEigenvalue sorting(m_ritz_val.data(), m_nev); - ind = sorting.index(); - } - break; - default: - throw std::invalid_argument("unsupported sorting rule"); - } - - ComplexVector new_ritz_val(m_ncv); - ComplexMatrix new_ritz_vec(m_ncv, m_nev); - BoolArray new_ritz_conv(m_nev); - - for(int i = 0; i < m_nev; i++) - { - new_ritz_val[i] = m_ritz_val[ind[i]]; - new_ritz_vec.col(i).noalias() = m_ritz_vec.col(ind[i]); - new_ritz_conv[i] = m_ritz_conv[ind[i]]; - } - - m_ritz_val.swap(new_ritz_val); - m_ritz_vec.swap(new_ritz_vec); - m_ritz_conv.swap(new_ritz_conv); - } + typedef Eigen::Index Index; public: /// @@ -496,7 +139,7 @@ /// 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 impelemnts all the public member functions + /// define their own that implements all the public member functions /// 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. @@ -506,224 +149,9 @@ /// 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, int nev, int ncv) : - m_op(op), - m_n(m_op->rows()), - m_nev(nev), - m_ncv(ncv > m_n ? m_n : ncv), - m_nmatop(0), - m_niter(0), - m_info(NOT_COMPUTED), - m_near_0(TypeTraits::min() * Scalar(10)), - m_eps(Eigen::NumTraits::epsilon()), - m_eps23(Eigen::numext::pow(m_eps, Scalar(2.0) / 3)) - { - if(nev < 1 || nev > m_n - 2) - throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 2, n is the size of matrix"); - - if(ncv < nev + 2 || ncv > m_n) - throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix"); - } - - /// - /// Virtual destructor - /// - virtual ~GenEigsSolver() {} - - /// - /// Initializes the solver by providing an initial residual vector. - /// - /// \param init_resid Pointer to the initial residual vector. - /// - /// **Spectra** (and also **ARPACK**) uses an iterative algorithm - /// to find eigenvalues. This function allows the user to provide the initial - /// residual vector. - /// - void init(const Scalar* init_resid) - { - // Reset all matrices/vectors to zero - m_fac_V.resize(m_n, m_ncv); - m_fac_H.resize(m_ncv, m_ncv); - m_fac_f.resize(m_n); - m_ritz_val.resize(m_ncv); - m_ritz_vec.resize(m_ncv, m_nev); - m_ritz_est.resize(m_ncv); - m_ritz_conv.resize(m_nev); - - m_fac_V.setZero(); - m_fac_H.setZero(); - m_fac_f.setZero(); - m_ritz_val.setZero(); - m_ritz_vec.setZero(); - m_ritz_est.setZero(); - m_ritz_conv.setZero(); - - // Set the initial vector - Vector v(m_n); - std::copy(init_resid, init_resid + m_n, v.data()); - const Scalar vnorm = v.norm(); - if(vnorm < m_near_0) - throw std::invalid_argument("initial residual vector cannot be zero"); - v /= vnorm; - - Vector w(m_n); - m_op->perform_op(v.data(), w.data()); - m_nmatop++; - - m_fac_H(0, 0) = v.dot(w); - m_fac_f.noalias() = w - v * m_fac_H(0, 0); - m_fac_V.col(0).noalias() = v; - - // In some cases f is zero in exact arithmetics, but due to rounding errors - // it may contain tiny fluctuations. When this happens, we force f to be zero - if(m_fac_f.cwiseAbs().maxCoeff() < m_eps) - m_fac_f.setZero(); - } - - /// - /// Initializes the solver by providing a random initial residual vector. - /// - /// This overloaded function generates a random initial residual vector - /// (with a fixed random seed) for the algorithm. Elements in the vector - /// follow independent Uniform(-0.5, 0.5) distribution. - /// - void init() - { - SimpleRandom rng(0); - Vector init_resid = rng.random_vec(m_n); - init(init_resid.data()); - } - - /// - /// Conducts the major computation procedure. - /// - /// \param maxit Maximum number of iterations allowed in the algorithm. - /// \param tol Precision parameter for the calculated eigenvalues. - /// \param sort_rule Rule to sort the eigenvalues and eigenvectors. - /// Supported values are - /// `Spectra::LARGEST_MAGN`, `Spectra::LARGEST_REAL`, - /// `Spectra::LARGEST_IMAG`, `Spectra::SMALLEST_MAGN`, - /// `Spectra::SMALLEST_REAL` and `Spectra::SMALLEST_IMAG`, - /// for example `LARGEST_MAGN` indicates that eigenvalues - /// with largest magnitude come first. - /// Note that this argument is only used to - /// **sort** the final result, and the **selection** rule - /// (e.g. selecting the largest or smallest eigenvalues in the - /// full spectrum) is specified by the template parameter - /// `SelectionRule` of GenEigsSolver. - /// - /// \return Number of converged eigenvalues. - /// - int compute(int maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_MAGN) - { - // The m-step Arnoldi factorization - factorize_from(1, m_ncv, m_fac_f); - retrieve_ritzpair(); - // Restarting - int i, nconv = 0, nev_adj; - for(i = 0; i < maxit; i++) - { - nconv = num_converged(tol); - if(nconv >= m_nev) - break; - - nev_adj = nev_adjusted(nconv); - restart(nev_adj); - } - // Sorting results - sort_ritzpair(sort_rule); - - m_niter += i + 1; - m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING; - - return std::min(m_nev, nconv); - } - - /// - /// Returns the status of the computation. - /// The full list of enumeration values can be found in \ref Enumerations. - /// - int info() const { return m_info; } - - /// - /// Returns the number of iterations used in the computation. - /// - int num_iterations() const { return m_niter; } - - /// - /// Returns the number of matrix operations used in the computation. - /// - int num_operations() const { return m_nmatop; } - - /// - /// Returns the converged eigenvalues. - /// - /// \return A complex-valued vector containing the eigenvalues. - /// Returned vector type will be `Eigen::Vector, ...>`, depending on - /// the template parameter `Scalar` defined. - /// - ComplexVector eigenvalues() const - { - const int nconv = m_ritz_conv.cast().sum(); - ComplexVector res(nconv); - - if(!nconv) - return res; - - int j = 0; - for(int i = 0; i < m_nev; i++) - { - if(m_ritz_conv[i]) - { - res[j] = m_ritz_val[i]; - j++; - } - } - - return res; - } - - /// - /// Returns the eigenvectors associated with the converged eigenvalues. - /// - /// \param nvec The number of eigenvectors to return. - /// - /// \return A complex-valued matrix containing the eigenvectors. - /// Returned matrix type will be `Eigen::Matrix, ...>`, - /// depending on the template parameter `Scalar` defined. - /// - ComplexMatrix eigenvectors(int nvec) const - { - const int nconv = m_ritz_conv.cast().sum(); - nvec = std::min(nvec, nconv); - ComplexMatrix res(m_n, nvec); - - if(!nvec) - return res; - - ComplexMatrix ritz_vec_conv(m_ncv, nvec); - int j = 0; - for(int i = 0; i < m_nev && j < nvec; i++) - { - if(m_ritz_conv[i]) - { - ritz_vec_conv.col(j).noalias() = m_ritz_vec.col(i); - j++; - } - } - - res.noalias() = m_fac_V * ritz_vec_conv; - - return res; - } - - /// - /// Returns all converged eigenvectors. - /// - ComplexMatrix eigenvectors() const - { - return eigenvectors(m_nev); - } + GenEigsSolver(OpType* op, Index nev, Index ncv) : + GenEigsBase(op, NULL, nev, ncv) + {} }; diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/Arnoldi.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/Arnoldi.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/Arnoldi.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/Arnoldi.h 2019-06-07 23:15:39.000000000 +0000 @@ -0,0 +1,283 @@ +// Copyright (C) 2018-2019 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 ARNOLDI_H +#define ARNOLDI_H + +#include +#include // std::sqrt +#include // std::invalid_argument +#include // std::stringstream + +#include "../MatOp/internal/ArnoldiOp.h" +#include "../Util/TypeTraits.h" +#include "../Util/SimpleRandom.h" +#include "UpperHessenbergQR.h" +#include "DoubleShiftQR.h" + +namespace Spectra { + + +// Arnoldi factorization A * V = V * H + f * e' +// A: n x n +// V: n x k +// H: k x k +// f: n x 1 +// e: [0, ..., 0, 1] +// V and H are allocated of dimension m, so the maximum value of k is m +template +class Arnoldi +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Matrix; + typedef Eigen::Matrix Vector; + typedef Eigen::Map MapMat; + typedef Eigen::Map MapVec; + typedef Eigen::Map MapConstMat; + typedef Eigen::Map MapConstVec; + +protected: + ArnoldiOpType m_op; // Operators for the Arnoldi factorization + + const Index m_n; // dimension of A + const Index m_m; // maximum dimension of subspace V + Index m_k; // current dimension of subspace V + + Matrix m_fac_V; // V matrix in the Arnoldi factorization + Matrix m_fac_H; // H matrix in the Arnoldi factorization + Vector m_fac_f; // residual in the Arnoldi factorization + Scalar m_beta; // ||f||, B-norm of f + + const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow + // ~= 1e-307 for the "double" type + const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type + + // Given orthonormal basis functions V, find a nonzero vector f such that V'Bf = 0 + // Assume that f has been properly allocated + void expand_basis(MapConstMat& V, const Index seed, Vector& f, Scalar& fnorm) + { + using std::sqrt; + + const Scalar thresh = m_eps * sqrt(Scalar(m_n)); + Vector Vf(V.cols()); + for(Index iter = 0; iter < 5; iter++) + { + // Randomly generate a new vector and orthogonalize it against V + SimpleRandom rng(seed + 123 * iter); + f.noalias() = rng.random_vec(m_n); + // f <- f - V * V'Bf, so that f is orthogonal to V in B-norm + m_op.trans_product(V, f, Vf); + f.noalias() -= V * Vf; + // fnorm <- ||f|| + fnorm = m_op.norm(f); + + // If fnorm is too close to zero, we try a new random vector, + // otherwise return the result + if(fnorm >= thresh) + return; + } + } + +public: + Arnoldi(const ArnoldiOpType& op, Index m) : + m_op(op), m_n(op.rows()), m_m(m), m_k(0), + m_near_0(TypeTraits::min() * Scalar(10)), + m_eps(Eigen::NumTraits::epsilon()) + {} + + virtual ~Arnoldi() {} + + // Const-reference to internal structures + const Matrix& matrix_V() const { return m_fac_V; } + const Matrix& matrix_H() const { return m_fac_H; } + const Vector& vector_f() const { return m_fac_f; } + Scalar f_norm() const { return m_beta; } + Index subspace_dim() const { return m_k; } + + // Initialize with an operator and an initial vector + void init(MapConstVec& v0, Index& op_counter) + { + m_fac_V.resize(m_n, m_m); + m_fac_H.resize(m_m, m_m); + m_fac_f.resize(m_n); + m_fac_H.setZero(); + + // Verify the initial vector + const Scalar v0norm = m_op.norm(v0); + if(v0norm < m_near_0) + throw std::invalid_argument("initial residual vector cannot be zero"); + + // Points to the first column of V + MapVec v(m_fac_V.data(), m_n); + + // Normalize + v.noalias() = v0 / v0norm; + + // Compute H and f + Vector w(m_n); + m_op.perform_op(v.data(), w.data()); + op_counter++; + + m_fac_H(0, 0) = m_op.inner_product(v, w); + m_fac_f.noalias() = w - v * m_fac_H(0, 0); + + // In some cases f is zero in exact arithmetics, but due to rounding errors + // it may contain tiny fluctuations. When this happens, we force f to be zero + if(m_fac_f.cwiseAbs().maxCoeff() < m_eps) + { + m_fac_f.setZero(); + m_beta = Scalar(0); + } else { + m_beta = m_op.norm(m_fac_f); + } + + // Indicate that this is a step-1 factorization + m_k = 1; + } + + // Arnoldi factorization starting from step-k + virtual void factorize_from(Index from_k, Index to_m, Index& op_counter) + { + using std::sqrt; + + if(to_m <= from_k) return; + + if(from_k > m_k) + { + std::stringstream msg; + msg << "Arnoldi: from_k (= " << from_k << + ") is larger than the current subspace dimension (= " << + m_k << ")"; + throw std::invalid_argument(msg.str()); + } + + const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n)); + + // Pre-allocate vectors + Vector Vf(to_m); + Vector w(m_n); + + // Keep the upperleft k x k submatrix of H and set other elements to 0 + m_fac_H.rightCols(m_m - from_k).setZero(); + m_fac_H.block(from_k, 0, m_m - from_k, from_k).setZero(); + + for(Index i = from_k; i <= to_m - 1; i++) + { + bool restart = false; + // If beta = 0, then the next V is not full rank + // We need to generate a new residual vector that is orthogonal + // to the current V, which we call a restart + if(m_beta < m_near_0) + { + MapConstMat V(m_fac_V.data(), m_n, i); // The first i columns + expand_basis(V, 2 * i, m_fac_f, m_beta); + restart = true; + } + + // v <- f / ||f|| + m_fac_V.col(i).noalias() = m_fac_f / m_beta; // The (i+1)-th column + + // Note that H[i+1, i] equals to the unrestarted beta + m_fac_H(i, i - 1) = restart ? Scalar(0) : m_beta; + + // w <- A * v, v = m_fac_V.col(i) + m_op.perform_op(&m_fac_V(0, i), w.data()); + op_counter++; + + const Index i1 = i + 1; + // First i+1 columns of V + MapConstMat Vs(m_fac_V.data(), m_n, i1); + // h = m_fac_H(0:i, i) + MapVec h(&m_fac_H(0, i), i1); + // h <- V'Bw + m_op.trans_product(Vs, w, h); + + // f <- w - V * h + m_fac_f.noalias() = w - Vs * h; + m_beta = m_op.norm(m_fac_f); + + if(m_beta > Scalar(0.717) * m_op.norm(h)) + continue; + + // f/||f|| is going to be the next column of V, so we need to test + // whether V'B(f/||f||) ~= 0 + m_op.trans_product(Vs, m_fac_f, Vf.head(i1)); + Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); + // If not, iteratively correct the residual + int count = 0; + while(count < 5 && ortho_err > m_eps * m_beta) + { + // There is an edge case: when beta=||f|| is close to zero, f mostly consists + // of noises of rounding errors, so the test [ortho_err < eps * beta] is very + // likely to fail. In particular, if beta=0, then the test is ensured to fail. + // Hence when this happens, we force f to be zero, and then restart in the + // next iteration. + if(m_beta < beta_thresh) + { + m_fac_f.setZero(); + m_beta = Scalar(0); + break; + } + + // f <- f - V * Vf + m_fac_f.noalias() -= Vs * Vf.head(i1); + // h <- h + Vf + h.noalias() += Vf.head(i1); + // beta <- ||f|| + m_beta = m_op.norm(m_fac_f); + + m_op.trans_product(Vs, m_fac_f, Vf.head(i1)); + ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); + count++; + } + } + + // Indicate that this is a step-m factorization + m_k = to_m; + } + + // Apply H -> Q'HQ, where Q is from a double shift QR decomposition + void compress_H(const DoubleShiftQR& decomp) + { + decomp.matrix_QtHQ(m_fac_H); + m_k -= 2; + } + + // Apply H -> Q'HQ, where Q is from an upper Hessenberg QR decomposition + void compress_H(const UpperHessenbergQR& decomp) + { + decomp.matrix_QtHQ(m_fac_H); + m_k--; + } + + // Apply V -> VQ and compute the new f. + // Should be called after compress_H(), since m_k is updated there. + // Only need to update the first k+1 columns of V + // The first (m - k + i) elements of the i-th column of Q are non-zero, + // and the rest are zero + void compress_V(const Matrix& Q) + { + Matrix Vs(m_n, m_k + 1); + for(Index i = 0; i < m_k; i++) + { + const Index nnz = m_m - m_k + i + 1; + MapConstVec q(&Q(0, i), nnz); + Vs.col(i).noalias() = m_fac_V.leftCols(nnz) * q; + } + Vs.col(m_k).noalias() = m_fac_V * Q.col(m_k); + m_fac_V.leftCols(m_k + 1).noalias() = Vs; + + Vector fk = m_fac_f * Q(m_m - 1, m_k - 1) + m_fac_V.col(m_k) * m_fac_H(m_k, m_k - 1); + m_fac_f.swap(fk); + m_beta = m_op.norm(m_fac_f); + } +}; + + +} // namespace Spectra + +#endif // ARNOLDI_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/BKLDLT.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/BKLDLT.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/BKLDLT.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/BKLDLT.h 2019-06-07 23:15:39.000000000 +0000 @@ -0,0 +1,501 @@ +// Copyright (C) 2019 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 BK_LDLT_H +#define BK_LDLT_H + +#include +#include +#include + +#include "../Util/CompInfo.h" + +namespace Spectra { + + +// Bunch-Kaufman LDLT decomposition +// References: +// 1. Bunch, J. R., & Kaufman, L. (1977). Some stable methods for calculating inertia and solving symmetric linear systems. +// Mathematics of computation, 31(137), 163-179. +// 2. Golub, G. H., & Van Loan, C. F. (2012). Matrix computations (Vol. 3). JHU press. Section 4.4. +// 3. Bunch-Parlett diagonal pivoting +template +class BKLDLT +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Matrix; + typedef Eigen::Matrix Vector; + typedef Eigen::Map MapVec; + typedef Eigen::Map MapConstVec; + + typedef Eigen::Matrix IntVector; + typedef Eigen::Ref GenericVector; + typedef Eigen::Ref GenericMatrix; + typedef const Eigen::Ref ConstGenericMatrix; + typedef const Eigen::Ref ConstGenericVector; + + Index m_n; + Vector m_data; // storage for a lower-triangular matrix + std::vector m_colptr; // pointers to columns + IntVector m_perm; // [-2, -1, 3, 1, 4, 5]: 0 <-> 2, 1 <-> 1, 2 <-> 3, 3 <-> 1, 4 <-> 4, 5 <-> 5 + std::vector< std::pair > m_permc; // compressed version of m_perm: [(0, 2), (2, 3), (3, 1)] + + bool m_computed; + int m_info; + + // Access to elements + // Pointer to the k-th column + Scalar* col_pointer(Index k) { return m_colptr[k]; } + // A[i, j] -> m_colptr[j][i - j], i >= j + Scalar& coeff(Index i, Index j) { return m_colptr[j][i - j]; } + const Scalar& coeff(Index i, Index j) const { return m_colptr[j][i - j]; } + // A[i, i] -> m_colptr[i][0] + Scalar& diag_coeff(Index i) { return m_colptr[i][0]; } + const Scalar& diag_coeff(Index i) const { return m_colptr[i][0]; } + + // Compute column pointers + void compute_pointer() + { + m_colptr.clear(); + m_colptr.reserve(m_n); + Scalar* head = m_data.data(); + + for(Index i = 0; i < m_n; i++) + { + m_colptr.push_back(head); + head += (m_n - i); + } + } + + // Copy mat - shift * I to m_data + void copy_data(ConstGenericMatrix& mat, int uplo, const Scalar& shift) + { + if(uplo == Eigen::Lower) + { + for(Index j = 0; j < m_n; j++) + { + const Scalar* begin = &mat.coeffRef(j, j); + const Index len = m_n - j; + std::copy(begin, begin + len, col_pointer(j)); + diag_coeff(j) -= shift; + } + } else { + Scalar* dest = m_data.data(); + for(Index i = 0; i < m_n; i++) + { + for(Index j = i; j < m_n; j++, dest++) + { + *dest = mat.coeff(i, j); + } + diag_coeff(i) -= shift; + } + } + } + + // Compute compressed permutations + void compress_permutation() + { + for(Index i = 0; i < m_n; i++) + { + // Recover the permutation action + const Index perm = (m_perm[i] >= 0) ? (m_perm[i]) : (-m_perm[i] - 1); + if(perm != i) + m_permc.push_back(std::make_pair(i, perm)); + } + } + + // Working on the A[k:end, k:end] submatrix + // Exchange k <-> r + // Assume r >= k + void pivoting_1x1(Index k, Index r) + { + // No permutation + if(k == r) + { + m_perm[k] = r; + return; + } + + // A[k, k] <-> A[r, r] + std::swap(diag_coeff(k), diag_coeff(r)); + + // A[(r+1):end, k] <-> A[(r+1):end, r] + std::swap_ranges(&coeff(r + 1, k), col_pointer(k + 1), &coeff(r + 1, r)); + + // A[(k+1):(r-1), k] <-> A[r, (k+1):(r-1)] + Scalar* src = &coeff(k + 1, k); + for(Index j = k + 1; j < r; j++, src++) + { + std::swap(*src, coeff(r, j)); + } + + m_perm[k] = r; + } + + // Working on the A[k:end, k:end] submatrix + // Exchange [k+1, k] <-> [r, p] + // Assume p >= k, r >= k+1 + void pivoting_2x2(Index k, Index r, Index p) + { + pivoting_1x1(k, p); + pivoting_1x1(k + 1, r); + + // A[k+1, k] <-> A[r, k] + std::swap(coeff(k + 1, k), coeff(r, k)); + + // Use negative signs to indicate a 2x2 block + // Also minus one to distinguish a negative zero from a positive zero + m_perm[k] = -m_perm[k] - 1; + m_perm[k + 1] = -m_perm[k + 1] - 1; + } + + // A[r1, c1:c2] <-> A[r2, c1:c2] + // Assume r2 >= r1 > c2 >= c1 + void interchange_rows(Index r1, Index r2, Index c1, Index c2) + { + if(r1 == r2) + return; + + for(Index j = c1; j <= c2; j++) + { + std::swap(coeff(r1, j), coeff(r2, j)); + } + } + + // lambda = |A[r, k]| = max{|A[k+1, k]|, ..., |A[end, k]|} + // Assume k < end + Scalar find_lambda(Index k, Index& r) + { + using std::abs; + + const Scalar* head = col_pointer(k); + const Scalar* end = col_pointer(k + 1); + Scalar lambda = abs(head[1]); + r = k + 1; + for(const Scalar* ptr = head + 2; ptr < end; ptr++) + { + const Scalar abs_elem = abs(*ptr); + if(lambda < abs_elem) + { + lambda = abs_elem; + r = k + (ptr - head); + } + } + + return lambda; + } + + // sigma = |A[p, r]| = max {|A[k, r]|, ..., |A[end, r]|} \ {A[r, r]} + // Assume k < r < end + Scalar find_sigma(Index k, Index r, Index& p) + { + using std::abs; + + // First search A[r+1, r], ..., A[end, r], which has the same task as find_lambda() + // If r == end, we skip this search + Scalar sigma = Scalar(-1); + if(r < m_n - 1) + sigma = find_lambda(r, p); + + // Then search A[k, r], ..., A[r-1, r], which maps to A[r, k], ..., A[r, r-1] + for(Index j = k; j < r; j++) + { + const Scalar abs_elem = abs(coeff(r, j)); + if(sigma < abs_elem) + { + sigma = abs_elem; + p = j; + } + } + + return sigma; + } + + // Generate permutations and apply to A + // Return true if the resulting pivoting is 1x1, and false if 2x2 + bool permutate_mat(Index k, const Scalar& alpha) + { + using std::abs; + + Index r = k, p = k; + const Scalar lambda = find_lambda(k, r); + + // If lambda=0, no need to interchange + if(lambda > Scalar(0)) + { + const Scalar abs_akk = abs(diag_coeff(k)); + // If |A[k, k]| >= alpha * lambda, no need to interchange + if(abs_akk < alpha * lambda) + { + const Scalar sigma = find_sigma(k, r, p); + + // If sigma * |A[k, k]| >= alpha * lambda^2, no need to interchange + if(sigma * abs_akk < alpha * lambda * lambda) + { + if(abs_akk >= alpha * sigma) + { + // Permutation on A + pivoting_1x1(k, r); + + // Permutation on L + interchange_rows(k, r, 0, k - 1); + return true; + } else { + // Permutation on A + // [r, p] and [p, r] are symmetric, but we need to make sure + // p >= k and r >= k+1, so it is safe to always make r > p + // One exception is when min{r,p} == k+1, in which case we make + // r = k+1, so that only one permutation needs to be performed + const Index rp_min = std::min(r, p); + const Index rp_max = std::max(r, p); + if(rp_min == k + 1) + { + r = rp_min; p = rp_max; + } else { + r = rp_max; p = rp_min; + } + pivoting_2x2(k, r, p); + // Permutation on L + interchange_rows(k, p, 0, k - 1); + interchange_rows(k + 1, r, 0, k - 1); + return false; + } + } + } + } + + return true; + } + + // E = [e11, e12] + // [e21, e22] + // Overwrite E with inv(E) + void inverse_inplace_2x2(Scalar& e11, Scalar& e21, Scalar& e22) const + { + // inv(E) = [d11, d12], d11 = e22/delta, d21 = -e21/delta, d22 = e11/delta + // [d21, d22] + const Scalar delta = e11 * e22 - e21 * e21; + std::swap(e11, e22); + e11 /= delta; + e22 /= delta; + e21 = -e21 / delta; + } + + // Return value is the status, SUCCESSFUL/NUMERICAL_ISSUE + int gaussian_elimination_1x1(Index k) + { + // D = 1 / A[k, k] + const Scalar akk = diag_coeff(k); + // Return NUMERICAL_ISSUE if not invertible + if(akk == Scalar(0)) + return NUMERICAL_ISSUE; + + diag_coeff(k) = 1.0 / akk; + + // B -= l * l' / A[k, k], B := A[(k+1):end, (k+1):end], l := L[(k+1):end, k] + Scalar* lptr = col_pointer(k) + 1; + const Index ldim = m_n - k - 1; + MapVec l(lptr, ldim); + for(Index j = 0; j < ldim; j++) + { + MapVec(col_pointer(j + k + 1), ldim - j).noalias() -= (lptr[j] / akk) * l.tail(ldim - j); + } + + // l /= A[k, k] + l /= akk; + + return SUCCESSFUL; + } + + // Return value is the status, SUCCESSFUL/NUMERICAL_ISSUE + int gaussian_elimination_2x2(Index k) + { + // D = inv(E) + Scalar& e11 = diag_coeff(k); + Scalar& e21 = coeff(k + 1, k); + Scalar& e22 = diag_coeff(k + 1); + // Return NUMERICAL_ISSUE if not invertible + if(e11 * e22 - e21 * e21 == Scalar(0)) + return NUMERICAL_ISSUE; + + inverse_inplace_2x2(e11, e21, e22); + + // X = l * inv(E), l := L[(k+2):end, k:(k+1)] + Scalar* l1ptr = &coeff(k + 2, k); + Scalar* l2ptr = &coeff(k + 2, k + 1); + const Index ldim = m_n - k - 2; + MapVec l1(l1ptr, ldim), l2(l2ptr, ldim); + + Eigen::Matrix X(ldim, 2); + X.col(0).noalias() = l1 * e11 + l2 * e21; + X.col(1).noalias() = l1 * e21 + l2 * e22; + + // B -= l * inv(E) * l' = X * l', B = A[(k+2):end, (k+2):end] + for(Index j = 0; j < ldim; j++) + { + MapVec(col_pointer(j + k + 2), ldim - j).noalias() -= (X.col(0).tail(ldim - j) * l1ptr[j] + X.col(1).tail(ldim - j) * l2ptr[j]); + } + + // l = X + l1.noalias() = X.col(0); + l2.noalias() = X.col(1); + + return SUCCESSFUL; + } + +public: + BKLDLT() : + m_n(0), m_computed(false), m_info(NOT_COMPUTED) + {} + + // Factorize mat - shift * I + BKLDLT(ConstGenericMatrix& mat, int uplo = Eigen::Lower, const Scalar& shift = Scalar(0)) : + m_n(mat.rows()), m_computed(false), m_info(NOT_COMPUTED) + { + compute(mat, uplo, shift); + } + + void compute(ConstGenericMatrix& mat, int uplo = Eigen::Lower, const Scalar& shift = Scalar(0)) + { + using std::abs; + + m_n = mat.rows(); + if(m_n != mat.cols()) + throw std::invalid_argument("BKLDLT: matrix must be square"); + + m_perm.setLinSpaced(m_n, 0, m_n - 1); + m_permc.clear(); + + // Copy data + m_data.resize((m_n * (m_n + 1)) / 2); + compute_pointer(); + copy_data(mat, uplo, shift); + + const Scalar alpha = (1.0 + std::sqrt(17.0)) / 8.0; + Index k = 0; + for(k = 0; k < m_n - 1; k++) + { + // 1. Interchange rows and columns of A, and save the result to m_perm + bool is_1x1 = permutate_mat(k, alpha); + + // 2. Gaussian elimination + if(is_1x1) + { + m_info = gaussian_elimination_1x1(k); + } else { + m_info = gaussian_elimination_2x2(k); + k++; + } + + // 3. Check status + if(m_info != SUCCESSFUL) + break; + } + // Invert the last 1x1 block if it exists + if(k == m_n - 1) + { + const Scalar akk = diag_coeff(k); + if(akk == Scalar(0)) + m_info = NUMERICAL_ISSUE; + + diag_coeff(k) = Scalar(1) / diag_coeff(k); + } + + compress_permutation(); + + m_computed = true; + } + + // Solve Ax=b + void solve_inplace(GenericVector b) const + { + if(!m_computed) + throw std::logic_error("BKLDLT: need to call compute() first"); + + // PAP' = LDL' + // 1. b -> Pb + Scalar* x = b.data(); + MapVec res(x, m_n); + Index npermc = m_permc.size(); + for(Index i = 0; i < npermc; i++) + { + std::swap(x[m_permc[i].first], x[m_permc[i].second]); + } + + // 2. Lz = Pb + // If m_perm[end] < 0, then end with m_n - 3, otherwise end with m_n - 2 + const Index end = (m_perm[m_n - 1] < 0) ? (m_n - 3) : (m_n - 2); + for(Index i = 0; i <= end; i++) + { + const Index b1size = m_n - i - 1; + const Index b2size = b1size - 1; + if(m_perm[i] >= 0) + { + MapConstVec l(&coeff(i + 1, i), b1size); + res.segment(i + 1, b1size).noalias() -= l * x[i]; + } else { + MapConstVec l1(&coeff(i + 2, i), b2size); + MapConstVec l2(&coeff(i + 2, i + 1), b2size); + res.segment(i + 2, b2size).noalias() -= (l1 * x[i] + l2 * x[i + 1]); + i++; + } + } + + // 3. Dw = z + for(Index i = 0; i < m_n; i++) + { + const Scalar e11 = diag_coeff(i); + if(m_perm[i] >= 0) + { + x[i] *= e11; + } else { + const Scalar e21 = coeff(i + 1, i), e22 = diag_coeff(i + 1); + const Scalar wi = x[i] * e11 + x[i + 1] * e21; + x[i + 1] = x[i] * e21 + x[i + 1] * e22; + x[i] = wi; + i++; + } + } + + // 4. L'y = w + // If m_perm[end] < 0, then start with m_n - 3, otherwise start with m_n - 2 + Index i = (m_perm[m_n - 1] < 0) ? (m_n - 3) : (m_n - 2); + for(; i >= 0; i--) + { + const Index ldim = m_n - i - 1; + MapConstVec l(&coeff(i + 1, i), ldim); + x[i] -= res.segment(i + 1, ldim).dot(l); + + if(m_perm[i] < 0) + { + MapConstVec l2(&coeff(i + 1, i - 1), ldim); + x[i - 1] -= res.segment(i + 1, ldim).dot(l2); + i--; + } + } + + // 5. x = P'y + for(Index i = npermc - 1; i >= 0; i--) + { + std::swap(x[m_permc[i].first], x[m_permc[i].second]); + } + } + + Vector solve(ConstGenericVector& b) const + { + Vector res = b; + solve_inplace(res); + return res; + } + + int info() const { return m_info; } +}; + + +} // namespace Spectra + +#endif // BK_LDLT_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/DoubleShiftQR.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/DoubleShiftQR.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/DoubleShiftQR.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/DoubleShiftQR.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -22,13 +22,12 @@ class DoubleShiftQR { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Matrix3X; typedef Eigen::Matrix Vector; typedef Eigen::Array IntArray; - typedef typename Matrix::Index Index; - typedef Eigen::Ref GenericMatrix; typedef const Eigen::Ref ConstGenericMatrix; @@ -160,7 +159,7 @@ // P = I - 2 * u * u' = P' // PX = X - 2 * u * (u'X) - void apply_PX(GenericMatrix X, Index stride, Index u_ind) + void apply_PX(GenericMatrix X, Index stride, Index u_ind) const { const Index nr = m_ref_nr.coeff(u_ind); if(nr == 1) @@ -198,7 +197,7 @@ // x is a pointer to a vector // Px = x - 2 * dot(x, u) * u - void apply_PX(Scalar* x, Index u_ind) + void apply_PX(Scalar* x, Index u_ind) const { const Index nr = m_ref_nr.coeff(u_ind); if(nr == 1) @@ -218,7 +217,7 @@ } // XP = X - 2 * (X * u) * u' - void apply_XP(GenericMatrix X, Index stride, Index u_ind) + void apply_XP(GenericMatrix X, Index stride, Index u_ind) const { const Index nr = m_ref_nr.coeff(u_ind); if(nr == 1) @@ -268,7 +267,7 @@ m_computed(false) {} - DoubleShiftQR(ConstGenericMatrix& mat, Scalar s, Scalar t) : + DoubleShiftQR(ConstGenericMatrix& mat, const Scalar& s, const Scalar& t) : m_n(mat.rows()), m_mat_H(m_n, m_n), m_shift_s(s), @@ -333,17 +332,17 @@ m_computed = true; } - const Matrix& matrix_QtHQ() const + void matrix_QtHQ(Matrix& dest) const { if(!m_computed) throw std::logic_error("DoubleShiftQR: need to call compute() first"); - return m_mat_H; + dest.noalias() = m_mat_H; } // Q = P0 * P1 * ... // Q'y = P_{n-2} * ... * P1 * P0 * y - void apply_QtY(Vector& y) + void apply_QtY(Vector& y) const { if(!m_computed) throw std::logic_error("DoubleShiftQR: need to call compute() first"); @@ -358,7 +357,7 @@ // Q = P0 * P1 * ... // YQ = Y * P0 * P1 * ... - void apply_YQ(GenericMatrix Y) + void apply_YQ(GenericMatrix Y) const { if(!m_computed) throw std::logic_error("DoubleShiftQR: need to call compute() first"); diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/Lanczos.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/Lanczos.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/Lanczos.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/Lanczos.h 2019-06-07 23:15:39.000000000 +0000 @@ -0,0 +1,170 @@ +// Copyright (C) 2018-2019 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 LANCZOS_H +#define LANCZOS_H + +#include +#include // std::sqrt +#include // std::invalid_argument +#include // std::stringstream + +#include "Arnoldi.h" + +namespace Spectra { + + +// Lanczos factorization A * V = V * H + f * e' +// A: n x n +// V: n x k +// H: k x k +// f: n x 1 +// e: [0, ..., 0, 1] +// V and H are allocated of dimension m, so the maximum value of k is m +template +class Lanczos: public Arnoldi +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Matrix; + typedef Eigen::Matrix Vector; + typedef Eigen::Map MapMat; + typedef Eigen::Map MapVec; + typedef Eigen::Map MapConstMat; + typedef Eigen::Map MapConstVec; + + using Arnoldi::m_op; + using Arnoldi::m_n; + using Arnoldi::m_m; + using Arnoldi::m_k; + using Arnoldi::m_fac_V; + using Arnoldi::m_fac_H; + using Arnoldi::m_fac_f; + using Arnoldi::m_beta; + using Arnoldi::m_near_0; + using Arnoldi::m_eps; + +public: + Lanczos(const ArnoldiOpType& op, Index m) : + Arnoldi(op, m) + {} + + // Lanczos factorization starting from step-k + void factorize_from(Index from_k, Index to_m, Index& op_counter) + { + using std::sqrt; + + if(to_m <= from_k) return; + + if(from_k > m_k) + { + std::stringstream msg; + msg << "Lanczos: from_k (= " << from_k << + ") is larger than the current subspace dimension (= " << + m_k << ")"; + throw std::invalid_argument(msg.str()); + } + + const Scalar beta_thresh = m_eps * sqrt(Scalar(m_n)); + + // Pre-allocate vectors + Vector Vf(to_m); + Vector w(m_n); + + // Keep the upperleft k x k submatrix of H and set other elements to 0 + m_fac_H.rightCols(m_m - from_k).setZero(); + m_fac_H.block(from_k, 0, m_m - from_k, from_k).setZero(); + + for(Index i = from_k; i <= to_m - 1; i++) + { + bool restart = false; + // If beta = 0, then the next V is not full rank + // We need to generate a new residual vector that is orthogonal + // to the current V, which we call a restart + if(m_beta < m_near_0) + { + MapConstMat V(m_fac_V.data(), m_n, i); // The first i columns + this->expand_basis(V, 2 * i, m_fac_f, m_beta); + restart = true; + } + + // v <- f / ||f|| + MapVec v(&m_fac_V(0, i), m_n); // The (i+1)-th column + v.noalias() = m_fac_f / m_beta; + + // Note that H[i+1, i] equals to the unrestarted beta + m_fac_H(i, i - 1) = restart ? Scalar(0) : m_beta; + + // w <- A * v + m_op.perform_op(v.data(), w.data()); + op_counter++; + + // H[i+1, i+1] = = v'Bw + m_fac_H(i - 1, i) = m_fac_H(i, i - 1); // Due to symmetry + m_fac_H(i, i) = m_op.inner_product(v, w); + + // f <- w - V * V'Bw = w - H[i+1, i] * V{i} - H[i+1, i+1] * V{i+1} + // If restarting, we know that H[i+1, i] = 0 + if(restart) + m_fac_f.noalias() = w - m_fac_H(i, i) * v; + else + m_fac_f.noalias() = w - m_fac_H(i, i - 1) * m_fac_V.col(i - 1) - m_fac_H(i, i) * v; + + m_beta = m_op.norm(m_fac_f); + + // f/||f|| is going to be the next column of V, so we need to test + // whether V'B(f/||f||) ~= 0 + const Index i1 = i + 1; + MapMat Vs(m_fac_V.data(), m_n, i1); // The first (i+1) columns + m_op.trans_product(Vs, m_fac_f, Vf.head(i1)); + Scalar ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); + // If not, iteratively correct the residual + int count = 0; + while(count < 5 && ortho_err > m_eps * m_beta) + { + // There is an edge case: when beta=||f|| is close to zero, f mostly consists + // of noises of rounding errors, so the test [ortho_err < eps * beta] is very + // likely to fail. In particular, if beta=0, then the test is ensured to fail. + // Hence when this happens, we force f to be zero, and then restart in the + // next iteration. + if(m_beta < beta_thresh) + { + m_fac_f.setZero(); + m_beta = Scalar(0); + break; + } + + // f <- f - V * Vf + m_fac_f.noalias() -= Vs * Vf.head(i1); + // h <- h + Vf + m_fac_H(i - 1, i) += Vf[i - 1]; + m_fac_H(i, i - 1) = m_fac_H(i - 1, i); + m_fac_H(i, i) += Vf[i]; + // beta <- ||f|| + m_beta = m_op.norm(m_fac_f); + + m_op.trans_product(Vs, m_fac_f, Vf.head(i1)); + ortho_err = Vf.head(i1).cwiseAbs().maxCoeff(); + count++; + } + } + + // Indicate that this is a step-m factorization + m_k = to_m; + } + + // Apply H -> Q'HQ, where Q is from a tridiagonal QR decomposition + void compress_H(const TridiagQR& decomp) + { + decomp.matrix_QtHQ(m_fac_H); + m_k--; + } +}; + + +} // namespace Spectra + +#endif // LANCZOS_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/TridiagEigen.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/TridiagEigen.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/TridiagEigen.h 2018-05-22 15:29:11.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/TridiagEigen.h 2019-06-07 23:15:39.000000000 +0000 @@ -2,7 +2,7 @@ // // Copyright (C) 2008-2010 Gael Guennebaud // Copyright (C) 2010 Jitse Niesen -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -24,14 +24,13 @@ class TridiagEigen { private: + typedef Eigen::Index Index; // For convenience in adapting the tridiagonal_qr_step() function typedef Scalar RealScalar; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef typename Matrix::Index Index; - typedef Eigen::Ref GenericMatrix; typedef const Eigen::Ref ConstGenericMatrix; @@ -152,7 +151,7 @@ Index end = m_n - 1; Index start = 0; - int iter = 0; // total number of iterations + Index iter = 0; // total number of iterations int info = 0; // 0 for success, 1 for failure const Scalar considerAsZero = TypeTraits::min(); diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/UpperHessenbergEigen.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/UpperHessenbergEigen.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/UpperHessenbergEigen.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/UpperHessenbergEigen.h 2019-06-07 23:15:39.000000000 +0000 @@ -2,7 +2,7 @@ // // Copyright (C) 2008 Gael Guennebaud // Copyright (C) 2010,2012 Jitse Niesen -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -22,11 +22,10 @@ class UpperHessenbergEigen { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef typename Matrix::Index Index; - typedef Eigen::Ref GenericMatrix; typedef const Eigen::Ref ConstGenericMatrix; diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/UpperHessenbergQR.h r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/UpperHessenbergQR.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/LinAlg/UpperHessenbergQR.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/LinAlg/UpperHessenbergQR.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -16,9 +16,21 @@ /// +/// \defgroup Internals Internal Classes +/// +/// Classes for internal use. May be useful to developers. +/// + +/// +/// \ingroup Internals +/// @{ +/// + +/// /// \defgroup LinearAlgebra Linear Algebra /// /// A number of classes for linear algebra operations. +/// /// /// \ingroup LinearAlgebra @@ -32,22 +44,23 @@ class UpperHessenbergQR { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; typedef Eigen::Matrix RowVector; typedef Eigen::Array Array; - typedef typename Matrix::Index Index; - typedef Eigen::Ref GenericMatrix; typedef const Eigen::Ref ConstGenericMatrix; + Matrix m_mat_T; + protected: Index m_n; - Matrix m_mat_T; // Gi = [ cos[i] sin[i]] // [-sin[i] cos[i]] // Q = G1 * G2 * ... * G_{n-1} + Scalar m_shift; Array m_rot_cos; Array m_rot_sin; bool m_computed; @@ -87,16 +100,21 @@ public: /// - /// Default constructor. Computation can + /// Constructor to preallocate memory. Computation can /// be performed later by calling the compute() method. /// - UpperHessenbergQR() : - m_n(0), m_computed(false) + UpperHessenbergQR(Index size) : + m_n(size), + m_rot_cos(m_n - 1), + m_rot_sin(m_n - 1), + m_computed(false) {} /// /// Constructor to create an object that performs and stores the - /// QR decomposition of an upper Hessenberg matrix `mat`. + /// QR decomposition of an upper Hessenberg matrix `mat`, with an + /// optional shift: \f$H-sI=QR\f$. Here \f$H\f$ stands for the matrix + /// `mat`, and \f$s\f$ is the shift. /// /// \param mat Matrix type can be `Eigen::Matrix` (e.g. /// `Eigen::MatrixXd` and `Eigen::MatrixXf`), or its mapped version @@ -104,23 +122,24 @@ /// Only the upper triangular and the lower subdiagonal parts of /// the matrix are used. /// - UpperHessenbergQR(ConstGenericMatrix& mat) : + UpperHessenbergQR(ConstGenericMatrix& mat, const Scalar& shift = Scalar(0)) : m_n(mat.rows()), - m_mat_T(m_n, m_n), + m_shift(shift), m_rot_cos(m_n - 1), m_rot_sin(m_n - 1), m_computed(false) { - compute(mat); + compute(mat, shift); } /// - /// We have virtual functions, so need a virtual destructor + /// Virtual destructor. /// - virtual ~UpperHessenbergQR(){}; + virtual ~UpperHessenbergQR() {}; /// - /// Conduct the QR factorization of an upper Hessenberg matrix. + /// Conduct the QR factorization of an upper Hessenberg matrix with + /// an optional shift. /// /// \param mat Matrix type can be `Eigen::Matrix` (e.g. /// `Eigen::MatrixXd` and `Eigen::MatrixXf`), or its mapped version @@ -128,18 +147,20 @@ /// Only the upper triangular and the lower subdiagonal parts of /// the matrix are used. /// - virtual void compute(ConstGenericMatrix& mat) + virtual void compute(ConstGenericMatrix& mat, const Scalar& shift = Scalar(0)) { m_n = mat.rows(); if(m_n != mat.cols()) throw std::invalid_argument("UpperHessenbergQR: matrix must be square"); + m_shift = shift; m_mat_T.resize(m_n, m_n); m_rot_cos.resize(m_n - 1); m_rot_sin.resize(m_n - 1); - // Make a copy of mat + // Make a copy of mat - s * I std::copy(mat.data(), mat.data() + mat.size(), m_mat_T.data()); + m_mat_T.diagonal().array() -= m_shift; Scalar xi, xj, r, c, s; Scalar *Tii, *ptr; @@ -202,13 +223,13 @@ } /// - /// Overwrite `dest` with the \f$RQ\f$ matrix, the multiplication of \f$R\f$ and \f$Q\f$. - /// \f$RQ\f$ is an upper Hessenberg matrix. + /// Overwrite `dest` with \f$Q'HQ = RQ + sI\f$, where \f$H\f$ is the input matrix `mat`, + /// and \f$s\f$ is the shift. The result is an upper Hessenberg matrix. /// /// \param mat The matrix to be overwritten, whose type should be `Eigen::Matrix`, /// depending on the template parameter `Scalar` defined. /// - virtual void matrix_RQ(Matrix& dest) const + virtual void matrix_QtHQ(Matrix& dest) const { if(!m_computed) throw std::logic_error("UpperHessenbergQR: need to call compute() first"); @@ -217,6 +238,7 @@ dest.resize(m_n, m_n); std::copy(m_mat_T.data(), m_mat_T.data() + m_mat_T.size(), dest.data()); + // Compute the RQ matrix const Index n1 = m_n - 1; for(Index i = 0; i < n1; i++) { @@ -240,6 +262,9 @@ dest.block(0, i, i + 2, 1) = c * Yi - s * dest.block(0, i + 1, i + 2, 1); dest.block(0, i + 1, i + 2, 1) = s * Yi + c * dest.block(0, i + 1, i + 2, 1); */ } + + // Add the shift to the diagonal + dest.diagonal().array() += m_shift; } /// @@ -437,6 +462,8 @@ /// +/// \ingroup LinearAlgebra +/// /// Perform the QR decomposition of a tridiagonal matrix, a special /// case of upper Hessenberg matrices. /// @@ -460,16 +487,18 @@ public: /// - /// Default constructor. Computation can + /// Constructor to preallocate memory. Computation can /// be performed later by calling the compute() method. /// - TridiagQR() : - UpperHessenbergQR() + TridiagQR(Index size) : + UpperHessenbergQR(size) {} /// /// Constructor to create an object that performs and stores the - /// QR decomposition of a tridiagonal matrix `mat`. + /// QR decomposition of an upper Hessenberg matrix `mat`, with an + /// optional shift: \f$H-sI=QR\f$. Here \f$H\f$ stands for the matrix + /// `mat`, and \f$s\f$ is the shift. /// /// \param mat Matrix type can be `Eigen::Matrix` (e.g. /// `Eigen::MatrixXd` and `Eigen::MatrixXf`), or its mapped version @@ -477,14 +506,15 @@ /// Only the major- and sub- diagonal parts of /// the matrix are used. /// - TridiagQR(ConstGenericMatrix& mat) : - UpperHessenbergQR() + TridiagQR(ConstGenericMatrix& mat, const Scalar& shift = Scalar(0)) : + UpperHessenbergQR(mat.rows()) { - this->compute(mat); + this->compute(mat, shift); } /// - /// Conduct the QR factorization of a tridiagonal matrix. + /// Conduct the QR factorization of a tridiagonal matrix with an + /// optional shift. /// /// \param mat Matrix type can be `Eigen::Matrix` (e.g. /// `Eigen::MatrixXd` and `Eigen::MatrixXf`), or its mapped version @@ -492,12 +522,13 @@ /// Only the major- and sub- diagonal parts of /// the matrix are used. /// - void compute(ConstGenericMatrix& mat) + void compute(ConstGenericMatrix& mat, const Scalar& shift = Scalar(0)) { this->m_n = mat.rows(); if(this->m_n != mat.cols()) throw std::invalid_argument("TridiagQR: matrix must be square"); + this->m_shift = shift; m_T_diag.resize(this->m_n); m_T_lsub.resize(this->m_n - 1); m_T_usub.resize(this->m_n - 1); @@ -505,7 +536,7 @@ this->m_rot_cos.resize(this->m_n - 1); this->m_rot_sin.resize(this->m_n - 1); - m_T_diag.noalias() = mat.diagonal(); + m_T_diag.array() = mat.diagonal().array() - this->m_shift; m_T_lsub.noalias() = mat.diagonal(-1); m_T_usub.noalias() = m_T_lsub; @@ -582,13 +613,13 @@ } /// - /// Overwrite `dest` with the \f$RQ\f$ matrix, the multiplication of \f$R\f$ and \f$Q\f$. - /// \f$RQ\f$ is an upper Hessenberg matrix. + /// Overwrite `dest` with \f$Q'HQ = RQ + sI\f$, where \f$H\f$ is the input matrix `mat`, + /// and \f$s\f$ is the shift. The result is a tridiagonal matrix. /// /// \param mat The matrix to be overwritten, whose type should be `Eigen::Matrix`, /// depending on the template parameter `Scalar` defined. /// - void matrix_RQ(Matrix& dest) const + void matrix_QtHQ(Matrix& dest) const { if(!this->m_computed) throw std::logic_error("TridiagQR: need to call compute() first"); @@ -600,6 +631,7 @@ // The upper diagonal refers to m_T_usub // The 2nd upper subdiagonal will be zero in RQ + // Compute the RQ matrix // [m11 m12] points to RQ[i:(i+1), i:(i+1)] // [0 m22] // @@ -622,9 +654,16 @@ // Copy the lower subdiagonal to upper subdiagonal dest.diagonal(1).noalias() = dest.diagonal(-1); + + // Add the shift to the diagonal + dest.diagonal().array() += this->m_shift; } }; +/// +/// @} +/// + } // namespace Spectra diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseCholesky.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseCholesky.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseCholesky.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseCholesky.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -10,6 +10,7 @@ #include #include #include +#include "../Util/CompInfo.h" namespace Spectra { @@ -26,15 +27,15 @@ class DenseCholesky { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; - typedef const Eigen::Ref ConstGenericMatrix; - const int m_n; + const Index m_n; Eigen::LLT m_decomp; int m_info; // status of the decomposition @@ -42,19 +43,18 @@ /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseCholesky(ConstGenericMatrix& mat_) : - m_n(mat_.rows()), - m_info(NOT_COMPUTED) + DenseCholesky(ConstGenericMatrix& mat) : + m_n(mat.rows()), m_info(NOT_COMPUTED) { - if(mat_.rows() != mat_.cols()) + if(mat.rows() != mat.cols()) throw std::invalid_argument("DenseCholesky: matrix must be square"); - m_decomp.compute(mat_); + m_decomp.compute(mat); m_info = (m_decomp.info() == Eigen::Success) ? SUCCESSFUL : NUMERICAL_ISSUE; @@ -63,11 +63,11 @@ /// /// Returns the number of rows of the underlying matrix. /// - int rows() const { return m_n; } + Index rows() const { return m_n; } /// /// Returns the number of columns of the underlying matrix. /// - int cols() const { return m_n; } + Index cols() const { return m_n; } /// /// Returns the status of the computation. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenComplexShiftSolve.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenComplexShiftSolve.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenComplexShiftSolve.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenComplexShiftSolve.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -26,11 +26,12 @@ class DenseGenComplexShiftSolve { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; + typedef const Eigen::Ref ConstGenericMatrix; typedef std::complex Complex; typedef Eigen::Matrix ComplexMatrix; @@ -38,10 +39,8 @@ typedef Eigen::PartialPivLU ComplexSolver; - typedef const Eigen::Ref ConstGenericMatrix; - - const MapConstMat m_mat; - const int m_n; + ConstGenericMatrix m_mat; + const Index m_n; ComplexSolver m_solver; ComplexVector m_x_cache; @@ -49,27 +48,26 @@ /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseGenComplexShiftSolve(ConstGenericMatrix& mat_) : - m_mat(mat_.data(), mat_.rows(), mat_.cols()), - m_n(mat_.rows()) + DenseGenComplexShiftSolve(ConstGenericMatrix& mat) : + m_mat(mat), m_n(mat.rows()) { - if(mat_.rows() != mat_.cols()) + if(mat.rows() != mat.cols()) throw std::invalid_argument("DenseGenComplexShiftSolve: matrix must be square"); } /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_n; } + Index rows() const { return m_n; } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_n; } + Index cols() const { return m_n; } /// /// Set the complex shift \f$\sigma\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenMatProd.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenMatProd.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenMatProd.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenMatProd.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -30,37 +30,36 @@ class DenseGenMatProd { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; - typedef const Eigen::Ref ConstGenericMatrix; - const MapConstMat m_mat; + ConstGenericMatrix m_mat; public: /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseGenMatProd(ConstGenericMatrix& mat_) : - m_mat(mat_.data(), mat_.rows(), mat_.cols()) + DenseGenMatProd(ConstGenericMatrix& mat) : + m_mat(mat) {} /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_mat.rows(); } + Index rows() const { return m_mat.rows(); } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_mat.cols(); } + Index cols() const { return m_mat.cols(); } /// /// Perform the matrix-vector multiplication operation \f$y=Ax\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenRealShiftSolve.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenRealShiftSolve.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseGenRealShiftSolve.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseGenRealShiftSolve.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -25,43 +25,41 @@ class DenseGenRealShiftSolve { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; - typedef const Eigen::Ref ConstGenericMatrix; - const MapConstMat m_mat; - const int m_n; + ConstGenericMatrix m_mat; + const Index m_n; Eigen::PartialPivLU m_solver; public: /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseGenRealShiftSolve(ConstGenericMatrix& mat_) : - m_mat(mat_.data(), mat_.rows(), mat_.cols()), - m_n(mat_.rows()) + DenseGenRealShiftSolve(ConstGenericMatrix& mat) : + m_mat(mat), m_n(mat.rows()) { - if(mat_.rows() != mat_.cols()) + if(mat.rows() != mat.cols()) throw std::invalid_argument("DenseGenRealShiftSolve: matrix must be square"); } /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_n; } + Index rows() const { return m_n; } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_n; } + Index cols() const { return m_n; } /// /// Set the real shift \f$\sigma\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseSymMatProd.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseSymMatProd.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseSymMatProd.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseSymMatProd.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -23,37 +23,36 @@ class DenseSymMatProd { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; - typedef const Eigen::Ref ConstGenericMatrix; - const MapConstMat m_mat; + ConstGenericMatrix m_mat; public: /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseSymMatProd(ConstGenericMatrix& mat_) : - m_mat(mat_.data(), mat_.rows(), mat_.cols()) + DenseSymMatProd(ConstGenericMatrix& mat) : + m_mat(mat) {} /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_mat.rows(); } + Index rows() const { return m_mat.rows(); } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_mat.cols(); } + Index cols() const { return m_mat.cols(); } /// /// Perform the matrix-vector multiplication operation \f$y=Ax\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseSymShiftSolve.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseSymShiftSolve.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/DenseSymShiftSolve.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/DenseSymShiftSolve.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -8,9 +8,11 @@ #define DENSE_SYM_SHIFT_SOLVE_H #include -#include #include +#include "../LinAlg/BKLDLT.h" +#include "../Util/CompInfo.h" + namespace Spectra { @@ -25,50 +27,50 @@ class DenseSymShiftSolve { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; - typedef Eigen::Map MapConstMat; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; - typedef const Eigen::Ref ConstGenericMatrix; - const MapConstMat m_mat; + ConstGenericMatrix m_mat; const int m_n; - Eigen::LDLT m_solver; + BKLDLT m_solver; public: /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** matrix object, whose type can be + /// \param mat An **Eigen** matrix object, whose type can be /// `Eigen::Matrix` (e.g. `Eigen::MatrixXd` and /// `Eigen::MatrixXf`), or its mapped version /// (e.g. `Eigen::Map`). /// - DenseSymShiftSolve(ConstGenericMatrix& mat_) : - m_mat(mat_.data(), mat_.rows(), mat_.cols()), - m_n(mat_.rows()) + DenseSymShiftSolve(ConstGenericMatrix& mat) : + m_mat(mat), m_n(mat.rows()) { - if(mat_.rows() != mat_.cols()) + if(mat.rows() != mat.cols()) throw std::invalid_argument("DenseSymShiftSolve: matrix must be square"); } /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_n; } + Index rows() const { return m_n; } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_n; } + Index cols() const { return m_n; } /// /// Set the real shift \f$\sigma\f$. /// void set_shift(Scalar sigma) { - m_solver.compute(m_mat - sigma * Matrix::Identity(m_n, m_n)); + m_solver.compute(m_mat, Uplo, sigma); + if(m_solver.info() != SUCCESSFUL) + throw std::invalid_argument("DenseSymShiftSolve: factorization failed with the given shift"); } /// diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/ArnoldiOp.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/ArnoldiOp.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/ArnoldiOp.h 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/ArnoldiOp.h 2019-06-07 23:15:39.000000000 +0000 @@ -0,0 +1,155 @@ +// Copyright (C) 2018-2019 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 ARNOLDI_OP_H +#define ARNOLDI_OP_H + +#include +#include // std::sqrt + +namespace Spectra { + + +/// +/// \ingroup Internals +/// @{ +/// + +/// +/// \defgroup Operators Operators +/// +/// Different types of operators. +/// + +/// +/// \ingroup Operators +/// +/// Operators used in the Arnoldi factorization. +/// +template +class ArnoldiOp +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Vector; + + OpType& m_op; + BOpType& m_Bop; + Vector m_cache; + +public: + ArnoldiOp(OpType* op, BOpType* Bop) : + m_op(*op), m_Bop(*Bop), m_cache(op->rows()) + {} + + inline Index rows() const { return m_op.rows(); } + + // In generalized eigenvalue problem Ax=lambda*Bx, define the inner product to be = x'By. + // For regular eigenvalue problems, it is the usual inner product = x'y + + // Compute = x'By + // x and y are two vectors + template + Scalar inner_product(const Arg1& x, const Arg2& y) + { + m_Bop.mat_prod(y.data(), m_cache.data()); + return x.dot(m_cache); + } + + // Compute res = = X'By + // X is a matrix, y is a vector, res is a vector + template + void trans_product(const Arg1& x, const Arg2& y, Eigen::Ref res) + { + m_Bop.mat_prod(y.data(), m_cache.data()); + res.noalias() = x.transpose() * m_cache; + } + + // B-norm of a vector, ||x||_B = sqrt(x'Bx) + template + Scalar norm(const Arg& x) + { + using std::sqrt; + return sqrt(inner_product(x, x)); + } + + // The "A" operator to generate the Krylov subspace + inline void perform_op(const Scalar* x_in, Scalar* y_out) + { + m_op.perform_op(x_in, y_out); + } +}; + + + +/// +/// \ingroup Operators +/// +/// Placeholder for the B-operator when \f$B = I\f$. +/// +class IdentityBOp {}; + + + +/// +/// \ingroup Operators +/// +/// Partial specialization for the case \f$B = I\f$. +/// +template +class ArnoldiOp +{ +private: + typedef Eigen::Index Index; + typedef Eigen::Matrix Vector; + + OpType& m_op; + +public: + ArnoldiOp(OpType* op, IdentityBOp* Bop) : + m_op(*op) + {} + + inline Index rows() const { return m_op.rows(); } + + // Compute = x'y + // x and y are two vectors + template + Scalar inner_product(const Arg1& x, const Arg2& y) const + { + return x.dot(y); + } + + // Compute res = = X'y + // X is a matrix, y is a vector, res is a vector + template + void trans_product(const Arg1& x, const Arg2& y, Eigen::Ref res) const + { + res.noalias() = x.transpose() * y; + } + + // B-norm of a vector. For regular eigenvalue problems it is simply the L2 norm + template + Scalar norm(const Arg& x) + { + return x.norm(); + } + + // The "A" operator to generate the Krylov subspace + inline void perform_op(const Scalar* x_in, Scalar* y_out) + { + m_op.perform_op(x_in, y_out); + } +}; + +/// +/// @} +/// + + +} // namespace Spectra + +#endif // ARNOLDI_OP_H diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/SymGEigsCholeskyOp.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/SymGEigsCholeskyOp.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/SymGEigsCholeskyOp.h 2018-05-22 13:09:59.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/SymGEigsCholeskyOp.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -15,7 +15,7 @@ /// -/// \ingroup MatOp +/// \ingroup Operators /// /// This class defines the matrix operation for generalized eigen solver in the /// Cholesky decomposition mode. It calculates \f$y=L^{-1}A(L')^{-1}x\f$ for any @@ -28,6 +28,7 @@ class SymGEigsCholeskyOp { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; @@ -49,11 +50,11 @@ /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_Bop.rows(); } + Index rows() const { return m_Bop.rows(); } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_Bop.rows(); } + Index cols() const { return m_Bop.rows(); } /// /// Perform the matrix operation \f$y=L^{-1}A(L')^{-1}x\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/SymGEigsRegInvOp.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/SymGEigsRegInvOp.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/internal/SymGEigsRegInvOp.h 2018-05-22 13:09:59.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/internal/SymGEigsRegInvOp.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2017-2018 Yixuan Qiu +// Copyright (C) 2017-2019 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 @@ -15,7 +15,7 @@ /// -/// \ingroup MatOp +/// \ingroup Operators /// /// This class defines the matrix operation for generalized eigen solver in the /// regular inverse mode. This class is intended for internal use. @@ -26,6 +26,7 @@ class SymGEigsRegInvOp { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Matrix; typedef Eigen::Matrix Vector; @@ -47,11 +48,11 @@ /// /// Return the number of rows of the underlying matrix. /// - int rows() const { return m_Bop.rows(); } + Index rows() const { return m_Bop.rows(); } /// /// Return the number of columns of the underlying matrix. /// - int cols() const { return m_Bop.rows(); } + Index cols() const { return m_Bop.rows(); } /// /// Perform the matrix operation \f$y=B^{-1}Ax\f$. diff -Nru r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/SparseCholesky.h r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/SparseCholesky.h --- r-cran-rspectra-0.13-1/inst/include/Spectra/MatOp/SparseCholesky.h 2018-03-04 01:06:19.000000000 +0000 +++ r-cran-rspectra-0.15-0/inst/include/Spectra/MatOp/SparseCholesky.h 2019-06-07 23:15:39.000000000 +0000 @@ -1,4 +1,4 @@ -// Copyright (C) 2016-2018 Yixuan Qiu +// Copyright (C) 2016-2019 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 @@ -11,7 +11,7 @@ #include #include #include -#include +#include "../Util/CompInfo.h" namespace Spectra { @@ -28,12 +28,14 @@ class SparseCholesky { private: + typedef Eigen::Index Index; typedef Eigen::Matrix Vector; typedef Eigen::Map MapConstVec; typedef Eigen::Map MapVec; typedef Eigen::SparseMatrix SparseMatrix; + typedef const Eigen::Ref ConstGenericSparseMatrix; - const int m_n; + const Index m_n; Eigen::SimplicialLLT m_decomp; int m_info; // status of the decomposition @@ -41,16 +43,17 @@ /// /// Constructor to create the matrix operation object. /// - /// \param mat_ An **Eigen** sparse matrix object, whose type is - /// `Eigen::SparseMatrix`. + /// \param mat An **Eigen** sparse matrix object, whose type can be + /// `Eigen::SparseMatrix

Large-Scale Eigenvalue Decomposition and SVD with RSpectra

Yixuan Qiu

2018-05-22

2019-06-10

Matrix in Sparse Format

Truncated SVD

Performance