src/Eigenvalues/ArpackSelfAdjointEigenSolver.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 David Harmon <dharmon@gmail.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
0011 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
0012
0013 #include "../../../../Eigen/Dense"
0014
0015 namespace Eigen {
0016
0017 namespace internal {
0018   template<typename Scalar, typename RealScalar> struct arpack_wrapper;
0019   template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> struct OP;
0020 }
0021
0022
0023
0024 template<typename MatrixType, typename MatrixSolver=SimplicialLLT<MatrixType>, bool BisSPD=false>
0025 class ArpackGeneralizedSelfAdjointEigenSolver
0026 {
0027 public:
0028   //typedef typename MatrixSolver::MatrixType MatrixType;
0029
0030   /** \brief Scalar type for matrices of type \p MatrixType. */
0031   typedef typename MatrixType::Scalar Scalar;
0032   typedef typename MatrixType::Index Index;
0033
0034   /** \brief Real scalar type for \p MatrixType.
0035    *
0036    * This is just \c Scalar if #Scalar is real (e.g., \c float or
0037    * \c Scalar), and the type of the real part of \c Scalar if #Scalar is
0038    * complex.
0039    */
0040   typedef typename NumTraits<Scalar>::Real RealScalar;
0041
0042   /** \brief Type for vector of eigenvalues as returned by eigenvalues().
0043    *
0044    * This is a column vector with entries of type #RealScalar.
0045    * The length of the vector is the size of \p nbrEigenvalues.
0046    */
0047   typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
0048
0049   /** \brief Default constructor.
0050    *
0051    * The default constructor is for cases in which the user intends to
0052    * perform decompositions via compute().
0053    *
0054    */
0055   ArpackGeneralizedSelfAdjointEigenSolver()
0056    : m_eivec(),
0057      m_eivalues(),
0058      m_isInitialized(false),
0059      m_eigenvectorsOk(false),
0060      m_nbrConverged(0),
0061      m_nbrIterations(0)
0062   { }
0063
0064   /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix.
0065    *
0066    * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
0067    *    computed. By default, the upper triangular part is used, but can be changed
0068    *    through the template parameter.
0069    * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem.
0070    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
0071    *    Must be less than the size of the input matrix, or an error is returned.
0072    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
0073    *    respective meanings to find the largest magnitude , smallest magnitude,
0074    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
0075    *    value can contain floating point value in string form, in which case the
0076    *    eigenvalues closest to this value will be found.
0077    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0078    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
0079    *    means machine precision.
0080    *
0081    * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar)
0082    * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if
0083    * \p options equals #ComputeEigenvectors.
0084    *
0085    */
0086   ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B,
0087                                           Index nbrEigenvalues, std::string eigs_sigma="LM",
0088                                int options=ComputeEigenvectors, RealScalar tol=0.0)
0089     : m_eivec(),
0090       m_eivalues(),
0091       m_isInitialized(false),
0092       m_eigenvectorsOk(false),
0093       m_nbrConverged(0),
0094       m_nbrIterations(0)
0095   {
0096     compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
0097   }
0098
0099   /** \brief Constructor; computes eigenvalues of given matrix.
0100    *
0101    * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
0102    *    computed. By default, the upper triangular part is used, but can be changed
0103    *    through the template parameter.
0104    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
0105    *    Must be less than the size of the input matrix, or an error is returned.
0106    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
0107    *    respective meanings to find the largest magnitude , smallest magnitude,
0108    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
0109    *    value can contain floating point value in string form, in which case the
0110    *    eigenvalues closest to this value will be found.
0111    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0112    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
0113    *    means machine precision.
0114    *
0115    * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar)
0116    * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if
0117    * \p options equals #ComputeEigenvectors.
0118    *
0119    */
0120
0121   ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A,
0122                                           Index nbrEigenvalues, std::string eigs_sigma="LM",
0123                                int options=ComputeEigenvectors, RealScalar tol=0.0)
0124     : m_eivec(),
0125       m_eivalues(),
0126       m_isInitialized(false),
0127       m_eigenvectorsOk(false),
0128       m_nbrConverged(0),
0129       m_nbrIterations(0)
0130   {
0131     compute(A, nbrEigenvalues, eigs_sigma, options, tol);
0132   }
0133
0134
0135   /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library.
0136    *
0137    * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
0138    * \param[in]  B  Selfadjoint matrix for generalized eigenvalues.
0139    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
0140    *    Must be less than the size of the input matrix, or an error is returned.
0141    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
0142    *    respective meanings to find the largest magnitude , smallest magnitude,
0143    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
0144    *    value can contain floating point value in string form, in which case the
0145    *    eigenvalues closest to this value will be found.
0146    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0147    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
0148    *    means machine precision.
0149    *
0150    * \returns    Reference to \c *this
0151    *
0152    * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK.  The eigenvalues()
0153    * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
0154    * then the eigenvectors are also computed and can be retrieved by
0155    * calling eigenvectors().
0156    *
0157    */
0158   ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, const MatrixType& B,
0159                                                    Index nbrEigenvalues, std::string eigs_sigma="LM",
0160                                         int options=ComputeEigenvectors, RealScalar tol=0.0);
0161
0162   /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library.
0163    *
0164    * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
0165    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
0166    *    Must be less than the size of the input matrix, or an error is returned.
0167    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
0168    *    respective meanings to find the largest magnitude , smallest magnitude,
0169    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
0170    *    value can contain floating point value in string form, in which case the
0171    *    eigenvalues closest to this value will be found.
0172    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0173    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
0174    *    means machine precision.
0175    *
0176    * \returns    Reference to \c *this
0177    *
0178    * This function computes the eigenvalues of \p A using ARPACK.  The eigenvalues()
0179    * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
0180    * then the eigenvectors are also computed and can be retrieved by
0181    * calling eigenvectors().
0182    *
0183    */
0184   ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A,
0185                                                    Index nbrEigenvalues, std::string eigs_sigma="LM",
0186                                         int options=ComputeEigenvectors, RealScalar tol=0.0);
0187
0188
0189   /** \brief Returns the eigenvectors of given matrix.
0190    *
0191    * \returns  A const reference to the matrix whose columns are the eigenvectors.
0192    *
0193    * \pre The eigenvectors have been computed before.
0194    *
0195    * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
0196    * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
0197    * eigenvectors are normalized to have (Euclidean) norm equal to one. If
0198    * this object was used to solve the eigenproblem for the selfadjoint
0199    * matrix \f$ A \f$, then the matrix returned by this function is the
0200    * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$.
0201    * For the generalized eigenproblem, the matrix returned is the solution \f$ A V = D B V \f$
0202    *
0203    * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp
0204    * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out
0205    *
0206    * \sa eigenvalues()
0207    */
0208   const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const
0209   {
0210     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
0211     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0212     return m_eivec;
0213   }
0214
0215   /** \brief Returns the eigenvalues of given matrix.
0216    *
0217    * \returns A const reference to the column vector containing the eigenvalues.
0218    *
0219    * \pre The eigenvalues have been computed before.
0220    *
0221    * The eigenvalues are repeated according to their algebraic multiplicity,
0222    * so there are as many eigenvalues as rows in the matrix. The eigenvalues
0223    * are sorted in increasing order.
0224    *
0225    * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
0226    * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
0227    *
0228    * \sa eigenvectors(), MatrixBase::eigenvalues()
0229    */
0230   const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
0231   {
0232     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
0233     return m_eivalues;
0234   }
0235
0236   /** \brief Computes the positive-definite square root of the matrix.
0237    *
0238    * \returns the positive-definite square root of the matrix
0239    *
0240    * \pre The eigenvalues and eigenvectors of a positive-definite matrix
0241    * have been computed before.
0242    *
0243    * The square root of a positive-definite matrix \f$ A \f$ is the
0244    * positive-definite matrix whose square equals \f$ A \f$. This function
0245    * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
0246    * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
0247    *
0248    * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp
0249    * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out
0250    *
0251    * \sa operatorInverseSqrt(),
0252    *     \ref MatrixFunctions_Module "MatrixFunctions Module"
0253    */
0254   Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const
0255   {
0256     eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0257     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0258     return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
0259   }
0260
0261   /** \brief Computes the inverse square root of the matrix.
0262    *
0263    * \returns the inverse positive-definite square root of the matrix
0264    *
0265    * \pre The eigenvalues and eigenvectors of a positive-definite matrix
0266    * have been computed before.
0267    *
0268    * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
0269    * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
0270    * cheaper than first computing the square root with operatorSqrt() and
0271    * then its inverse with MatrixBase::inverse().
0272    *
0273    * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp
0274    * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out
0275    *
0276    * \sa operatorSqrt(), MatrixBase::inverse(),
0277    *     \ref MatrixFunctions_Module "MatrixFunctions Module"
0278    */
0279   Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const
0280   {
0281     eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0282     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0283     return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
0284   }
0285
0286   /** \brief Reports whether previous computation was successful.
0287    *
0288    * \returns \c Success if computation was successful, \c NoConvergence otherwise.
0289    */
0290   ComputationInfo info() const
0291   {
0292     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
0293     return m_info;
0294   }
0295
0296   size_t getNbrConvergedEigenValues() const
0297   { return m_nbrConverged; }
0298
0299   size_t getNbrIterations() const
0300   { return m_nbrIterations; }
0301
0302 protected:
0303   Matrix<Scalar, Dynamic, Dynamic> m_eivec;
0304   Matrix<Scalar, Dynamic, 1> m_eivalues;
0305   ComputationInfo m_info;
0306   bool m_isInitialized;
0307   bool m_eigenvectorsOk;
0308
0309   size_t m_nbrConverged;
0310   size_t m_nbrIterations;
0311 };
0312
0313
0314
0315
0316
0317 template<typename MatrixType, typename MatrixSolver, bool BisSPD>
0318 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
0319     ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
0320 ::compute(const MatrixType& A, Index nbrEigenvalues,
0321           std::string eigs_sigma, int options, RealScalar tol)
0322 {
0323     MatrixType B(0,0);
0324     compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
0325
0326     return *this;
0327 }
0328
0329
0330 template<typename MatrixType, typename MatrixSolver, bool BisSPD>
0331 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
0332     ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
0333 ::compute(const MatrixType& A, const MatrixType& B, Index nbrEigenvalues,
0334           std::string eigs_sigma, int options, RealScalar tol)
0335 {
0336   eigen_assert(A.cols() == A.rows());
0337   eigen_assert(B.cols() == B.rows());
0338   eigen_assert(B.rows() == 0 || A.cols() == B.rows());
0339   eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0
0340             && (options & EigVecMask) != EigVecMask
0341             && "invalid option parameter");
0342
0343   bool isBempty = (B.rows() == 0) || (B.cols() == 0);
0344
0345   // For clarity, all parameters match their ARPACK name
0346   //
0347   // Always 0 on the first call
0348   //
0349   int ido = 0;
0350
0351   int n = (int)A.cols();
0352
0353   // User options: "LA", "SA", "SM", "LM", "BE"
0354   //
0355   char whch[3] = "LM";
0356
0357   // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
0358   //
0359   RealScalar sigma = 0.0;
0360
0361   if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
0362   {
0363       eigs_sigma[0] = toupper(eigs_sigma[0]);
0364       eigs_sigma[1] = toupper(eigs_sigma[1]);
0365
0366       // In the following special case we're going to invert the problem, since solving
0367       // for larger magnitude is much much faster
0368       // i.e., if 'SM' is specified, we're going to really use 'LM', the default
0369       //
0370       if (eigs_sigma.substr(0,2) != "SM")
0371       {
0372           whch[0] = eigs_sigma[0];
0373           whch[1] = eigs_sigma[1];
0374       }
0375   }
0376   else
0377   {
0378       eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
0379
0380       // If it's not scalar values, then the user may be explicitly
0381       // specifying the sigma value to cluster the evs around
0382       //
0383       sigma = atof(eigs_sigma.c_str());
0384
0385       // If atof fails, it returns 0.0, which is a fine default
0386       //
0387   }
0388
0389   // "I" means normal eigenvalue problem, "G" means generalized
0390   //
0391   char bmat[2] = "I";
0392   if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
0393       bmat[0] = 'G';
0394
0395   // Now we determine the mode to use
0396   //
0397   int mode = (bmat[0] == 'G') + 1;
0398   if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
0399   {
0400       // We're going to use shift-and-invert mode, and basically find
0401       // the largest eigenvalues of the inverse operator
0402       //
0403       mode = 3;
0404   }
0405
0406   // The user-specified number of eigenvalues/vectors to compute
0407   //
0408   int nev = (int)nbrEigenvalues;
0409
0410   // Allocate space for ARPACK to store the residual
0411   //
0412   Scalar *resid = new Scalar[n];
0413
0414   // Number of Lanczos vectors, must satisfy nev < ncv <= n
0415   // Note that this indicates that nev != n, and we cannot compute
0416   // all eigenvalues of a mtrix
0417   //
0418   int ncv = std::min(std::max(2*nev, 20), n);
0419
0420   // The working n x ncv matrix, also store the final eigenvectors (if computed)
0421   //
0422   Scalar *v = new Scalar[n*ncv];
0423   int ldv = n;
0424
0425   // Working space
0426   //
0427   Scalar *workd = new Scalar[3*n];
0428   int lworkl = ncv*ncv+8*ncv; // Must be at least this length
0429   Scalar *workl = new Scalar[lworkl];
0430
0431   int *iparam= new int[11];
0432   iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
0433   iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1)));
0434   iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
0435
0436   // Used during reverse communicate to notify where arrays start
0437   //
0438   int *ipntr = new int[11];
0439
0440   // Error codes are returned in here, initial value of 0 indicates a random initial
0441   // residual vector is used, any other values means resid contains the initial residual
0442   // vector, possibly from a previous run
0443   //
0444   int info = 0;
0445
0446   Scalar scale = 1.0;
0447   //if (!isBempty)
0448   //{
0449   //Scalar scale = B.norm() / std::sqrt(n);
0450   //scale = std::pow(2, std::floor(std::log(scale+1)));
0451   ////M /= scale;
0452   //for (size_t i=0; i<(size_t)B.outerSize(); i++)
0453   //    for (typename MatrixType::InnerIterator it(B, i); it; ++it)
0454   //        it.valueRef() /= scale;
0455   //}
0456
0457   MatrixSolver OP;
0458   if (mode == 1 || mode == 2)
0459   {
0460       if (!isBempty)
0461           OP.compute(B);
0462   }
0463   else if (mode == 3)
0464   {
0465       if (sigma == 0.0)
0466       {
0467           OP.compute(A);
0468       }
0469       else
0470       {
0471           // Note: We will never enter here because sigma must be 0.0
0472           //
0473           if (isBempty)
0474           {
0475             MatrixType AminusSigmaB(A);
0476             for (Index i=0; i<A.rows(); ++i)
0477                 AminusSigmaB.coeffRef(i,i) -= sigma;
0478
0479             OP.compute(AminusSigmaB);
0480           }
0481           else
0482           {
0483               MatrixType AminusSigmaB = A - sigma * B;
0484               OP.compute(AminusSigmaB);
0485           }
0486       }
0487   }
0488
0489   if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success)
0490       std::cout << "Error factoring matrix" << std::endl;
0491
0492   do
0493   {
0494     internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid,
0495                                                         &ncv, v, &ldv, iparam, ipntr, workd, workl,
0496                                                         &lworkl, &info);
0497
0498     if (ido == -1 || ido == 1)
0499     {
0500       Scalar *in  = workd + ipntr[0] - 1;
0501       Scalar *out = workd + ipntr[1] - 1;
0502
0503       if (ido == 1 && mode != 2)
0504       {
0505           Scalar *out2 = workd + ipntr[2] - 1;
0506           if (isBempty || mode == 1)
0507             Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
0508           else
0509             Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
0510
0511           in = workd + ipntr[2] - 1;
0512       }
0513
0514       if (mode == 1)
0515       {
0516         if (isBempty)
0517         {
0518           // OP = A
0519           //
0520           Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
0521         }
0522         else
0523         {
0524           // OP = L^{-1}AL^{-T}
0525           //
0526           internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
0527         }
0528       }
0529       else if (mode == 2)
0530       {
0531         if (ido == 1)
0532           Matrix<Scalar, Dynamic, 1>::Map(in, n)  = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
0533
0534         // OP = B^{-1} A
0535         //
0536         Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
0537       }
0538       else if (mode == 3)
0539       {
0540         // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
0541         // The B * in is already computed and stored at in if ido == 1
0542         //
0543         if (ido == 1 || isBempty)
0544           Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
0545         else
0546           Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
0547       }
0548     }
0549     else if (ido == 2)
0550     {
0551       Scalar *in  = workd + ipntr[0] - 1;
0552       Scalar *out = workd + ipntr[1] - 1;
0553
0554       if (isBempty || mode == 1)
0555         Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
0556       else
0557         Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
0558     }
0559   } while (ido != 99);
0560
0561   if (info == 1)
0562     m_info = NoConvergence;
0563   else if (info == 3)
0564     m_info = NumericalIssue;
0565   else if (info < 0)
0566     m_info = InvalidInput;
0567   else if (info != 0)
0568     eigen_assert(false && "Unknown ARPACK return value!");
0569   else
0570   {
0571     // Do we compute eigenvectors or not?
0572     //
0573     int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors;
0574
0575     // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
0576     //
0577     char howmny[2] = "A";
0578
0579     // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
0580     //
0581     int *select = new int[ncv];
0582
0583     // Final eigenvalues
0584     //
0585     m_eivalues.resize(nev, 1);
0586
0587     internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
0588                                                         &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
0589                                                         v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info);
0590
0591     if (info == -14)
0592       m_info = NoConvergence;
0593     else if (info != 0)
0594       m_info = InvalidInput;
0595     else
0596     {
0597       if (rvec)
0598       {
0599         m_eivec.resize(A.rows(), nev);
0600         for (int i=0; i<nev; i++)
0601           for (int j=0; j<n; j++)
0602             m_eivec(j,i) = v[i*n+j] / scale;
0603
0604         if (mode == 1 && !isBempty && BisSPD)
0605           internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
0606
0607         m_eigenvectorsOk = true;
0608       }
0609
0610       m_nbrIterations = iparam[2];
0611       m_nbrConverged  = iparam[4];
0612
0613       m_info = Success;
0614     }
0615
0616     delete[] select;
0617   }
0618
0619   delete[] v;
0620   delete[] iparam;
0621   delete[] ipntr;
0622   delete[] workd;
0623   delete[] workl;
0624   delete[] resid;
0625
0626   m_isInitialized = true;
0627
0628   return *this;
0629 }
0630
0631
0632 // Single precision
0633 //
0634 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which,
0635     int *nev, float *tol, float *resid, int *ncv,
0636     float *v, int *ldv, int *iparam, int *ipntr,
0637     float *workd, float *workl, int *lworkl,
0638     int *info);
0639
0640 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d,
0641     float *z, int *ldz, float *sigma,
0642     char *bmat, int *n, char *which, int *nev,
0643     float *tol, float *resid, int *ncv, float *v,
0644     int *ldv, int *iparam, int *ipntr, float *workd,
0645     float *workl, int *lworkl, int *ierr);
0646
0647 // Double precision
0648 //
0649 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
0650     int *nev, double *tol, double *resid, int *ncv,
0651     double *v, int *ldv, int *iparam, int *ipntr,
0652     double *workd, double *workl, int *lworkl,
0653     int *info);
0654
0655 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
0656     double *z, int *ldz, double *sigma,
0657     char *bmat, int *n, char *which, int *nev,
0658     double *tol, double *resid, int *ncv, double *v,
0659     int *ldv, int *iparam, int *ipntr, double *workd,
0660     double *workl, int *lworkl, int *ierr);
0661
0662
0663 namespace internal {
0664
0665 template<typename Scalar, typename RealScalar> struct arpack_wrapper
0666 {
0667   static inline void saupd(int *ido, char *bmat, int *n, char *which,
0668       int *nev, RealScalar *tol, Scalar *resid, int *ncv,
0669       Scalar *v, int *ldv, int *iparam, int *ipntr,
0670       Scalar *workd, Scalar *workl, int *lworkl, int *info)
0671   {
0672     EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
0673   }
0674
0675   static inline void seupd(int *rvec, char *All, int *select, Scalar *d,
0676       Scalar *z, int *ldz, RealScalar *sigma,
0677       char *bmat, int *n, char *which, int *nev,
0678       RealScalar *tol, Scalar *resid, int *ncv, Scalar *v,
0679       int *ldv, int *iparam, int *ipntr, Scalar *workd,
0680       Scalar *workl, int *lworkl, int *ierr)
0681   {
0682     EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
0683   }
0684 };
0685
0686 template <> struct arpack_wrapper<float, float>
0687 {
0688   static inline void saupd(int *ido, char *bmat, int *n, char *which,
0689       int *nev, float *tol, float *resid, int *ncv,
0690       float *v, int *ldv, int *iparam, int *ipntr,
0691       float *workd, float *workl, int *lworkl, int *info)
0692   {
0693     ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
0694   }
0695
0696   static inline void seupd(int *rvec, char *All, int *select, float *d,
0697       float *z, int *ldz, float *sigma,
0698       char *bmat, int *n, char *which, int *nev,
0699       float *tol, float *resid, int *ncv, float *v,
0700       int *ldv, int *iparam, int *ipntr, float *workd,
0701       float *workl, int *lworkl, int *ierr)
0702   {
0703     sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
0704         workd, workl, lworkl, ierr);
0705   }
0706 };
0707
0708 template <> struct arpack_wrapper<double, double>
0709 {
0710   static inline void saupd(int *ido, char *bmat, int *n, char *which,
0711       int *nev, double *tol, double *resid, int *ncv,
0712       double *v, int *ldv, int *iparam, int *ipntr,
0713       double *workd, double *workl, int *lworkl, int *info)
0714   {
0715     dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
0716   }
0717
0718   static inline void seupd(int *rvec, char *All, int *select, double *d,
0719       double *z, int *ldz, double *sigma,
0720       char *bmat, int *n, char *which, int *nev,
0721       double *tol, double *resid, int *ncv, double *v,
0722       int *ldv, int *iparam, int *ipntr, double *workd,
0723       double *workl, int *lworkl, int *ierr)
0724   {
0725     dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
0726         workd, workl, lworkl, ierr);
0727   }
0728 };
0729
0730
0731 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD>
0732 struct OP
0733 {
0734     static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out);
0735     static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs);
0736 };
0737
0738 template<typename MatrixSolver, typename MatrixType, typename Scalar>
0739 struct OP<MatrixSolver, MatrixType, Scalar, true>
0740 {
0741   static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
0742 {
0743     // OP = L^{-1} A L^{-T}  (B = LL^T)
0744     //
0745     // First solve L^T out = in
0746     //
0747     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
0748     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
0749
0750     // Then compute out = A out
0751     //
0752     Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n);
0753
0754     // Then solve L out = out
0755     //
0756     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
0757     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n));
0758 }
0759
0760   static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
0761 {
0762     // Solve L^T out = in
0763     //
0764     Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k));
0765     Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k);
0766 }
0767
0768 };
0769
0770 template<typename MatrixSolver, typename MatrixType, typename Scalar>
0771 struct OP<MatrixSolver, MatrixType, Scalar, false>
0772 {
0773   static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
0774 {
0775     eigen_assert(false && "Should never be in here...");
0776 }
0777
0778   static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
0779 {
0780     eigen_assert(false && "Should never be in here...");
0781 }
0782
0783 };
0784
0785 } // end namespace internal
0786
0787 } // end namespace Eigen
0788
0789 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
0790