Home | History | Annotate | Download | only in Eigenvalues
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 David Harmon <dharmon (at) gmail.com>
      5 //
      6 // Eigen is free software; you can redistribute it and/or
      7 // modify it under the terms of the GNU Lesser General Public
      8 // License as published by the Free Software Foundation; either
      9 // version 3 of the License, or (at your option) any later version.
     10 //
     11 // Alternatively, you can redistribute it and/or
     12 // modify it under the terms of the GNU General Public License as
     13 // published by the Free Software Foundation; either version 2 of
     14 // the License, or (at your option) any later version.
     15 //
     16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
     17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
     18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
     19 // GNU General Public License for more details.
     20 //
     21 // You should have received a copy of the GNU Lesser General Public
     22 // License and a copy of the GNU General Public License along with
     23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
     24 
     25 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
     26 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
     27 
     28 #include <Eigen/Dense>
     29 
     30 namespace Eigen {
     31 
     32 namespace internal {
     33   template<typename Scalar, typename RealScalar> struct arpack_wrapper;
     34   template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> struct OP;
     35 }
     36 
     37 
     38 
     39 template<typename MatrixType, typename MatrixSolver=SimplicialLLT<MatrixType>, bool BisSPD=false>
     40 class ArpackGeneralizedSelfAdjointEigenSolver
     41 {
     42 public:
     43   //typedef typename MatrixSolver::MatrixType MatrixType;
     44 
     45   /** \brief Scalar type for matrices of type \p MatrixType. */
     46   typedef typename MatrixType::Scalar Scalar;
     47   typedef typename MatrixType::Index Index;
     48 
     49   /** \brief Real scalar type for \p MatrixType.
     50    *
     51    * This is just \c Scalar if #Scalar is real (e.g., \c float or
     52    * \c Scalar), and the type of the real part of \c Scalar if #Scalar is
     53    * complex.
     54    */
     55   typedef typename NumTraits<Scalar>::Real RealScalar;
     56 
     57   /** \brief Type for vector of eigenvalues as returned by eigenvalues().
     58    *
     59    * This is a column vector with entries of type #RealScalar.
     60    * The length of the vector is the size of \p nbrEigenvalues.
     61    */
     62   typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
     63 
     64   /** \brief Default constructor.
     65    *
     66    * The default constructor is for cases in which the user intends to
     67    * perform decompositions via compute().
     68    *
     69    */
     70   ArpackGeneralizedSelfAdjointEigenSolver()
     71    : m_eivec(),
     72      m_eivalues(),
     73      m_isInitialized(false),
     74      m_eigenvectorsOk(false),
     75      m_nbrConverged(0),
     76      m_nbrIterations(0)
     77   { }
     78 
     79   /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix.
     80    *
     81    * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
     82    *    computed. By default, the upper triangular part is used, but can be changed
     83    *    through the template parameter.
     84    * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem.
     85    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
     86    *    Must be less than the size of the input matrix, or an error is returned.
     87    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
     88    *    respective meanings to find the largest magnitude , smallest magnitude,
     89    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
     90    *    value can contain floating point value in string form, in which case the
     91    *    eigenvalues closest to this value will be found.
     92    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
     93    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
     94    *    means machine precision.
     95    *
     96    * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar)
     97    * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if
     98    * \p options equals #ComputeEigenvectors.
     99    *
    100    */
    101   ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B,
    102                                           Index nbrEigenvalues, std::string eigs_sigma="LM",
    103                                int options=ComputeEigenvectors, RealScalar tol=0.0)
    104     : m_eivec(),
    105       m_eivalues(),
    106       m_isInitialized(false),
    107       m_eigenvectorsOk(false),
    108       m_nbrConverged(0),
    109       m_nbrIterations(0)
    110   {
    111     compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
    112   }
    113 
    114   /** \brief Constructor; computes eigenvalues of given matrix.
    115    *
    116    * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
    117    *    computed. By default, the upper triangular part is used, but can be changed
    118    *    through the template parameter.
    119    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
    120    *    Must be less than the size of the input matrix, or an error is returned.
    121    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
    122    *    respective meanings to find the largest magnitude , smallest magnitude,
    123    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
    124    *    value can contain floating point value in string form, in which case the
    125    *    eigenvalues closest to this value will be found.
    126    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
    127    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
    128    *    means machine precision.
    129    *
    130    * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar)
    131    * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if
    132    * \p options equals #ComputeEigenvectors.
    133    *
    134    */
    135 
    136   ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A,
    137                                           Index nbrEigenvalues, std::string eigs_sigma="LM",
    138                                int options=ComputeEigenvectors, RealScalar tol=0.0)
    139     : m_eivec(),
    140       m_eivalues(),
    141       m_isInitialized(false),
    142       m_eigenvectorsOk(false),
    143       m_nbrConverged(0),
    144       m_nbrIterations(0)
    145   {
    146     compute(A, nbrEigenvalues, eigs_sigma, options, tol);
    147   }
    148 
    149 
    150   /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library.
    151    *
    152    * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
    153    * \param[in]  B  Selfadjoint matrix for generalized eigenvalues.
    154    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
    155    *    Must be less than the size of the input matrix, or an error is returned.
    156    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
    157    *    respective meanings to find the largest magnitude , smallest magnitude,
    158    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
    159    *    value can contain floating point value in string form, in which case the
    160    *    eigenvalues closest to this value will be found.
    161    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
    162    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
    163    *    means machine precision.
    164    *
    165    * \returns    Reference to \c *this
    166    *
    167    * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK.  The eigenvalues()
    168    * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
    169    * then the eigenvectors are also computed and can be retrieved by
    170    * calling eigenvectors().
    171    *
    172    */
    173   ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, const MatrixType& B,
    174                                                    Index nbrEigenvalues, std::string eigs_sigma="LM",
    175                                         int options=ComputeEigenvectors, RealScalar tol=0.0);
    176 
    177   /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library.
    178    *
    179    * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
    180    * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
    181    *    Must be less than the size of the input matrix, or an error is returned.
    182    * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
    183    *    respective meanings to find the largest magnitude , smallest magnitude,
    184    *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
    185    *    value can contain floating point value in string form, in which case the
    186    *    eigenvalues closest to this value will be found.
    187    * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
    188    * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
    189    *    means machine precision.
    190    *
    191    * \returns    Reference to \c *this
    192    *
    193    * This function computes the eigenvalues of \p A using ARPACK.  The eigenvalues()
    194    * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
    195    * then the eigenvectors are also computed and can be retrieved by
    196    * calling eigenvectors().
    197    *
    198    */
    199   ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A,
    200                                                    Index nbrEigenvalues, std::string eigs_sigma="LM",
    201                                         int options=ComputeEigenvectors, RealScalar tol=0.0);
    202 
    203 
    204   /** \brief Returns the eigenvectors of given matrix.
    205    *
    206    * \returns  A const reference to the matrix whose columns are the eigenvectors.
    207    *
    208    * \pre The eigenvectors have been computed before.
    209    *
    210    * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
    211    * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
    212    * eigenvectors are normalized to have (Euclidean) norm equal to one. If
    213    * this object was used to solve the eigenproblem for the selfadjoint
    214    * matrix \f$ A \f$, then the matrix returned by this function is the
    215    * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$.
    216    * For the generalized eigenproblem, the matrix returned is the solution \f$ A V = D B V \f$
    217    *
    218    * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp
    219    * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out
    220    *
    221    * \sa eigenvalues()
    222    */
    223   const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const
    224   {
    225     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    226     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    227     return m_eivec;
    228   }
    229 
    230   /** \brief Returns the eigenvalues of given matrix.
    231    *
    232    * \returns A const reference to the column vector containing the eigenvalues.
    233    *
    234    * \pre The eigenvalues have been computed before.
    235    *
    236    * The eigenvalues are repeated according to their algebraic multiplicity,
    237    * so there are as many eigenvalues as rows in the matrix. The eigenvalues
    238    * are sorted in increasing order.
    239    *
    240    * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
    241    * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
    242    *
    243    * \sa eigenvectors(), MatrixBase::eigenvalues()
    244    */
    245   const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
    246   {
    247     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    248     return m_eivalues;
    249   }
    250 
    251   /** \brief Computes the positive-definite square root of the matrix.
    252    *
    253    * \returns the positive-definite square root of the matrix
    254    *
    255    * \pre The eigenvalues and eigenvectors of a positive-definite matrix
    256    * have been computed before.
    257    *
    258    * The square root of a positive-definite matrix \f$ A \f$ is the
    259    * positive-definite matrix whose square equals \f$ A \f$. This function
    260    * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
    261    * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
    262    *
    263    * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp
    264    * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out
    265    *
    266    * \sa operatorInverseSqrt(),
    267    *     \ref MatrixFunctions_Module "MatrixFunctions Module"
    268    */
    269   Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const
    270   {
    271     eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
    272     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    273     return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
    274   }
    275 
    276   /** \brief Computes the inverse square root of the matrix.
    277    *
    278    * \returns the inverse positive-definite square root of the matrix
    279    *
    280    * \pre The eigenvalues and eigenvectors of a positive-definite matrix
    281    * have been computed before.
    282    *
    283    * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
    284    * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
    285    * cheaper than first computing the square root with operatorSqrt() and
    286    * then its inverse with MatrixBase::inverse().
    287    *
    288    * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp
    289    * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out
    290    *
    291    * \sa operatorSqrt(), MatrixBase::inverse(),
    292    *     \ref MatrixFunctions_Module "MatrixFunctions Module"
    293    */
    294   Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const
    295   {
    296     eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
    297     eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
    298     return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
    299   }
    300 
    301   /** \brief Reports whether previous computation was successful.
    302    *
    303    * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
    304    */
    305   ComputationInfo info() const
    306   {
    307     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
    308     return m_info;
    309   }
    310 
    311   size_t getNbrConvergedEigenValues() const
    312   { return m_nbrConverged; }
    313 
    314   size_t getNbrIterations() const
    315   { return m_nbrIterations; }
    316 
    317 protected:
    318   Matrix<Scalar, Dynamic, Dynamic> m_eivec;
    319   Matrix<Scalar, Dynamic, 1> m_eivalues;
    320   ComputationInfo m_info;
    321   bool m_isInitialized;
    322   bool m_eigenvectorsOk;
    323 
    324   size_t m_nbrConverged;
    325   size_t m_nbrIterations;
    326 };
    327 
    328 
    329 
    330 
    331 
    332 template<typename MatrixType, typename MatrixSolver, bool BisSPD>
    333 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
    334     ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
    335 ::compute(const MatrixType& A, Index nbrEigenvalues,
    336           std::string eigs_sigma, int options, RealScalar tol)
    337 {
    338     MatrixType B(0,0);
    339     compute(A, B, nbrEigenvalues, eigs_sigma, options, tol);
    340 
    341     return *this;
    342 }
    343 
    344 
    345 template<typename MatrixType, typename MatrixSolver, bool BisSPD>
    346 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>&
    347     ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>
    348 ::compute(const MatrixType& A, const MatrixType& B, Index nbrEigenvalues,
    349           std::string eigs_sigma, int options, RealScalar tol)
    350 {
    351   eigen_assert(A.cols() == A.rows());
    352   eigen_assert(B.cols() == B.rows());
    353   eigen_assert(B.rows() == 0 || A.cols() == B.rows());
    354   eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0
    355             && (options & EigVecMask) != EigVecMask
    356             && "invalid option parameter");
    357 
    358   bool isBempty = (B.rows() == 0) || (B.cols() == 0);
    359 
    360   // For clarity, all parameters match their ARPACK name
    361   //
    362   // Always 0 on the first call
    363   //
    364   int ido = 0;
    365 
    366   int n = (int)A.cols();
    367 
    368   // User options: "LA", "SA", "SM", "LM", "BE"
    369   //
    370   char whch[3] = "LM";
    371 
    372   // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
    373   //
    374   RealScalar sigma = 0.0;
    375 
    376   if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
    377   {
    378       eigs_sigma[0] = toupper(eigs_sigma[0]);
    379       eigs_sigma[1] = toupper(eigs_sigma[1]);
    380 
    381       // In the following special case we're going to invert the problem, since solving
    382       // for larger magnitude is much much faster
    383       // i.e., if 'SM' is specified, we're going to really use 'LM', the default
    384       //
    385       if (eigs_sigma.substr(0,2) != "SM")
    386       {
    387           whch[0] = eigs_sigma[0];
    388           whch[1] = eigs_sigma[1];
    389       }
    390   }
    391   else
    392   {
    393       eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
    394 
    395       // If it's not scalar values, then the user may be explicitly
    396       // specifying the sigma value to cluster the evs around
    397       //
    398       sigma = atof(eigs_sigma.c_str());
    399 
    400       // If atof fails, it returns 0.0, which is a fine default
    401       //
    402   }
    403 
    404   // "I" means normal eigenvalue problem, "G" means generalized
    405   //
    406   char bmat[2] = "I";
    407   if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
    408       bmat[0] = 'G';
    409 
    410   // Now we determine the mode to use
    411   //
    412   int mode = (bmat[0] == 'G') + 1;
    413   if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
    414   {
    415       // We're going to use shift-and-invert mode, and basically find
    416       // the largest eigenvalues of the inverse operator
    417       //
    418       mode = 3;
    419   }
    420 
    421   // The user-specified number of eigenvalues/vectors to compute
    422   //
    423   int nev = (int)nbrEigenvalues;
    424 
    425   // Allocate space for ARPACK to store the residual
    426   //
    427   Scalar *resid = new Scalar[n];
    428 
    429   // Number of Lanczos vectors, must satisfy nev < ncv <= n
    430   // Note that this indicates that nev != n, and we cannot compute
    431   // all eigenvalues of a mtrix
    432   //
    433   int ncv = std::min(std::max(2*nev, 20), n);
    434 
    435   // The working n x ncv matrix, also store the final eigenvectors (if computed)
    436   //
    437   Scalar *v = new Scalar[n*ncv];
    438   int ldv = n;
    439 
    440   // Working space
    441   //
    442   Scalar *workd = new Scalar[3*n];
    443   int lworkl = ncv*ncv+8*ncv; // Must be at least this length
    444   Scalar *workl = new Scalar[lworkl];
    445 
    446   int *iparam= new int[11];
    447   iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
    448   iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1)));
    449   iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
    450 
    451   // Used during reverse communicate to notify where arrays start
    452   //
    453   int *ipntr = new int[11];
    454 
    455   // Error codes are returned in here, initial value of 0 indicates a random initial
    456   // residual vector is used, any other values means resid contains the initial residual
    457   // vector, possibly from a previous run
    458   //
    459   int info = 0;
    460 
    461   Scalar scale = 1.0;
    462   //if (!isBempty)
    463   //{
    464   //Scalar scale = B.norm() / std::sqrt(n);
    465   //scale = std::pow(2, std::floor(std::log(scale+1)));
    466   ////M /= scale;
    467   //for (size_t i=0; i<(size_t)B.outerSize(); i++)
    468   //    for (typename MatrixType::InnerIterator it(B, i); it; ++it)
    469   //        it.valueRef() /= scale;
    470   //}
    471 
    472   MatrixSolver OP;
    473   if (mode == 1 || mode == 2)
    474   {
    475       if (!isBempty)
    476           OP.compute(B);
    477   }
    478   else if (mode == 3)
    479   {
    480       if (sigma == 0.0)
    481       {
    482           OP.compute(A);
    483       }
    484       else
    485       {
    486           // Note: We will never enter here because sigma must be 0.0
    487           //
    488           if (isBempty)
    489           {
    490             MatrixType AminusSigmaB(A);
    491             for (Index i=0; i<A.rows(); ++i)
    492                 AminusSigmaB.coeffRef(i,i) -= sigma;
    493 
    494             OP.compute(AminusSigmaB);
    495           }
    496           else
    497           {
    498               MatrixType AminusSigmaB = A - sigma * B;
    499               OP.compute(AminusSigmaB);
    500           }
    501       }
    502   }
    503 
    504   if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success)
    505       std::cout << "Error factoring matrix" << std::endl;
    506 
    507   do
    508   {
    509     internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid,
    510                                                         &ncv, v, &ldv, iparam, ipntr, workd, workl,
    511                                                         &lworkl, &info);
    512 
    513     if (ido == -1 || ido == 1)
    514     {
    515       Scalar *in  = workd + ipntr[0] - 1;
    516       Scalar *out = workd + ipntr[1] - 1;
    517 
    518       if (ido == 1 && mode != 2)
    519       {
    520           Scalar *out2 = workd + ipntr[2] - 1;
    521           if (isBempty || mode == 1)
    522             Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
    523           else
    524             Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
    525 
    526           in = workd + ipntr[2] - 1;
    527       }
    528 
    529       if (mode == 1)
    530       {
    531         if (isBempty)
    532         {
    533           // OP = A
    534           //
    535           Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
    536         }
    537         else
    538         {
    539           // OP = L^{-1}AL^{-T}
    540           //
    541           internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
    542         }
    543       }
    544       else if (mode == 2)
    545       {
    546         if (ido == 1)
    547           Matrix<Scalar, Dynamic, 1>::Map(in, n)  = A * Matrix<Scalar, Dynamic, 1>::Map(in, n);
    548 
    549         // OP = B^{-1} A
    550         //
    551         Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
    552       }
    553       else if (mode == 3)
    554       {
    555         // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
    556         // The B * in is already computed and stored at in if ido == 1
    557         //
    558         if (ido == 1 || isBempty)
    559           Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
    560         else
    561           Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
    562       }
    563     }
    564     else if (ido == 2)
    565     {
    566       Scalar *in  = workd + ipntr[0] - 1;
    567       Scalar *out = workd + ipntr[1] - 1;
    568 
    569       if (isBempty || mode == 1)
    570         Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
    571       else
    572         Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
    573     }
    574   } while (ido != 99);
    575 
    576   if (info == 1)
    577     m_info = NoConvergence;
    578   else if (info == 3)
    579     m_info = NumericalIssue;
    580   else if (info < 0)
    581     m_info = InvalidInput;
    582   else if (info != 0)
    583     eigen_assert(false && "Unknown ARPACK return value!");
    584   else
    585   {
    586     // Do we compute eigenvectors or not?
    587     //
    588     int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors;
    589 
    590     // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
    591     //
    592     char howmny[2] = "A";
    593 
    594     // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
    595     //
    596     int *select = new int[ncv];
    597 
    598     // Final eigenvalues
    599     //
    600     m_eivalues.resize(nev, 1);
    601 
    602     internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
    603                                                         &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
    604                                                         v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info);
    605 
    606     if (info == -14)
    607       m_info = NoConvergence;
    608     else if (info != 0)
    609       m_info = InvalidInput;
    610     else
    611     {
    612       if (rvec)
    613       {
    614         m_eivec.resize(A.rows(), nev);
    615         for (int i=0; i<nev; i++)
    616           for (int j=0; j<n; j++)
    617             m_eivec(j,i) = v[i*n+j] / scale;
    618 
    619         if (mode == 1 && !isBempty && BisSPD)
    620           internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
    621 
    622         m_eigenvectorsOk = true;
    623       }
    624 
    625       m_nbrIterations = iparam[2];
    626       m_nbrConverged  = iparam[4];
    627 
    628       m_info = Success;
    629     }
    630 
    631     delete select;
    632   }
    633 
    634   delete v;
    635   delete iparam;
    636   delete ipntr;
    637   delete workd;
    638   delete workl;
    639   delete resid;
    640 
    641   m_isInitialized = true;
    642 
    643   return *this;
    644 }
    645 
    646 
    647 // Single precision
    648 //
    649 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which,
    650     int *nev, float *tol, float *resid, int *ncv,
    651     float *v, int *ldv, int *iparam, int *ipntr,
    652     float *workd, float *workl, int *lworkl,
    653     int *info);
    654 
    655 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d,
    656     float *z, int *ldz, float *sigma,
    657     char *bmat, int *n, char *which, int *nev,
    658     float *tol, float *resid, int *ncv, float *v,
    659     int *ldv, int *iparam, int *ipntr, float *workd,
    660     float *workl, int *lworkl, int *ierr);
    661 
    662 // Double precision
    663 //
    664 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
    665     int *nev, double *tol, double *resid, int *ncv,
    666     double *v, int *ldv, int *iparam, int *ipntr,
    667     double *workd, double *workl, int *lworkl,
    668     int *info);
    669 
    670 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
    671     double *z, int *ldz, double *sigma,
    672     char *bmat, int *n, char *which, int *nev,
    673     double *tol, double *resid, int *ncv, double *v,
    674     int *ldv, int *iparam, int *ipntr, double *workd,
    675     double *workl, int *lworkl, int *ierr);
    676 
    677 
    678 namespace internal {
    679 
    680 template<typename Scalar, typename RealScalar> struct arpack_wrapper
    681 {
    682   static inline void saupd(int *ido, char *bmat, int *n, char *which,
    683       int *nev, RealScalar *tol, Scalar *resid, int *ncv,
    684       Scalar *v, int *ldv, int *iparam, int *ipntr,
    685       Scalar *workd, Scalar *workl, int *lworkl, int *info)
    686   {
    687     EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
    688   }
    689 
    690   static inline void seupd(int *rvec, char *All, int *select, Scalar *d,
    691       Scalar *z, int *ldz, RealScalar *sigma,
    692       char *bmat, int *n, char *which, int *nev,
    693       RealScalar *tol, Scalar *resid, int *ncv, Scalar *v,
    694       int *ldv, int *iparam, int *ipntr, Scalar *workd,
    695       Scalar *workl, int *lworkl, int *ierr)
    696   {
    697     EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
    698   }
    699 };
    700 
    701 template <> struct arpack_wrapper<float, float>
    702 {
    703   static inline void saupd(int *ido, char *bmat, int *n, char *which,
    704       int *nev, float *tol, float *resid, int *ncv,
    705       float *v, int *ldv, int *iparam, int *ipntr,
    706       float *workd, float *workl, int *lworkl, int *info)
    707   {
    708     ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
    709   }
    710 
    711   static inline void seupd(int *rvec, char *All, int *select, float *d,
    712       float *z, int *ldz, float *sigma,
    713       char *bmat, int *n, char *which, int *nev,
    714       float *tol, float *resid, int *ncv, float *v,
    715       int *ldv, int *iparam, int *ipntr, float *workd,
    716       float *workl, int *lworkl, int *ierr)
    717   {
    718     sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
    719         workd, workl, lworkl, ierr);
    720   }
    721 };
    722 
    723 template <> struct arpack_wrapper<double, double>
    724 {
    725   static inline void saupd(int *ido, char *bmat, int *n, char *which,
    726       int *nev, double *tol, double *resid, int *ncv,
    727       double *v, int *ldv, int *iparam, int *ipntr,
    728       double *workd, double *workl, int *lworkl, int *info)
    729   {
    730     dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
    731   }
    732 
    733   static inline void seupd(int *rvec, char *All, int *select, double *d,
    734       double *z, int *ldz, double *sigma,
    735       char *bmat, int *n, char *which, int *nev,
    736       double *tol, double *resid, int *ncv, double *v,
    737       int *ldv, int *iparam, int *ipntr, double *workd,
    738       double *workl, int *lworkl, int *ierr)
    739   {
    740     dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
    741         workd, workl, lworkl, ierr);
    742   }
    743 };
    744 
    745 
    746 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD>
    747 struct OP
    748 {
    749     static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out);
    750     static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs);
    751 };
    752 
    753 template<typename MatrixSolver, typename MatrixType, typename Scalar>
    754 struct OP<MatrixSolver, MatrixType, Scalar, true>
    755 {
    756   static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
    757 {
    758     // OP = L^{-1} A L^{-T}  (B = LL^T)
    759     //
    760     // First solve L^T out = in
    761     //
    762     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n));
    763     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
    764 
    765     // Then compute out = A out
    766     //
    767     Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n);
    768 
    769     // Then solve L out = out
    770     //
    771     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n);
    772     Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n));
    773 }
    774 
    775   static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
    776 {
    777     // Solve L^T out = in
    778     //
    779     Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k));
    780     Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k);
    781 }
    782 
    783 };
    784 
    785 template<typename MatrixSolver, typename MatrixType, typename Scalar>
    786 struct OP<MatrixSolver, MatrixType, Scalar, false>
    787 {
    788   static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out)
    789 {
    790     eigen_assert(false && "Should never be in here...");
    791 }
    792 
    793   static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs)
    794 {
    795     eigen_assert(false && "Should never be in here...");
    796 }
    797 
    798 };
    799 
    800 } // end namespace internal
    801 
    802 } // end namespace Eigen
    803 
    804 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
    805 
    806