LCOV - code coverage report
Current view: top level - utils - arpack_wrapper.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 121 131 92.4 %
Date: 2013-05-24 Functions: 48 48 100.0 %
Branches: 629 1705 36.9 %

           Branch data     Line data    Source code
       1                 :            : /* This software is distributed under BSD 3-clause license (see LICENSE file).
       2                 :            :  *
       3                 :            :  * This code is based on Eigen3 ARPACK wrapper
       4                 :            :  * Copyright (c) 2012 David Harmon
       5                 :            :  *
       6                 :            :  * The code is relicensed under BSD 3-clause with the permission of its
       7                 :            :  * original author. Some changes were made afterwards by Sergey Lisitsyn.
       8                 :            :  */
       9                 :            : 
      10                 :            : #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
      11                 :            : #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
      12                 :            : 
      13                 :            : /* Tapkee includes */
      14                 :            : #include <tapkee/defines.hpp>
      15                 :            : /* End of Tapkee includes */
      16                 :            : 
      17                 :            : using Eigen::Matrix;
      18                 :            : using Eigen::Dynamic;
      19                 :            : using Eigen::ComputationInfo;
      20                 :            : using Eigen::ComputeEigenvectors;
      21                 :            : using Eigen::NoConvergence;
      22                 :            : using Eigen::NumericalIssue;
      23                 :            : #ifndef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
      24                 :            : using Eigen::InvalidInput;
      25                 :            : #endif
      26                 :            : using Eigen::Success;
      27                 :            : using Eigen::NumTraits;
      28                 :            : 
      29                 :            : //! Namespace that contains ARPACK routines
      30                 :            : namespace arpack_internal 
      31                 :            : {
      32                 :            :         template<typename Scalar, typename RealScalar> struct arpack_wrapper;
      33                 :            : }
      34                 :            : 
      35                 :            : template<class LMatrixType, class RMatrixType, class MatrixOperation, bool BisSPD=false>
      36 [ +  - ][ +  - ]:         27 : class ArpackGeneralizedSelfAdjointEigenSolver
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
      37                 :            : {
      38                 :            : public:
      39                 :            :   //typedef typename MatrixSolver::MatrixType MatrixType;
      40                 :            : 
      41                 :            :         /** \brief Scalar type for matrices of type \p MatrixType. */
      42                 :            :         typedef typename LMatrixType::Scalar Scalar;
      43                 :            :         typedef typename LMatrixType::Index Index;
      44                 :            : 
      45                 :            :         /** \brief Real scalar type for \p MatrixType.
      46                 :            :         *
      47                 :            :         * This is just \c Scalar if #Scalar is real (e.g., \c float or
      48                 :            :         * \c Scalar), and the type of the real part of \c Scalar if #Scalar is
      49                 :            :         * complex.
      50                 :            :         */
      51                 :            :         typedef typename NumTraits<Scalar>::Real RealScalar;
      52                 :            : 
      53                 :            :         /** \brief Type for vector of eigenvalues as returned by eigenvalues().
      54                 :            :         *
      55                 :            :         * This is a column vector with entries of type #RealScalar.
      56                 :            :         * The length of the vector is the size of \p nbrEigenvalues.
      57                 :            :         */
      58                 :            :         typedef typename Eigen::internal::plain_col_type<LMatrixType, RealScalar>::type RealVectorType;
      59                 :            : 
      60                 :            :         /** \brief Default constructor.
      61                 :            :         *
      62                 :            :         * The default constructor is for cases in which the user intends to
      63                 :            :         * perform decompositions via compute().
      64                 :            :         *
      65                 :            :         */
      66                 :            :         ArpackGeneralizedSelfAdjointEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false),
      67                 :            :                 m_eigenvectorsOk(false), m_nbrConverged(0), m_nbrIterations(0)
      68                 :            :         { 
      69                 :            :         }
      70                 :            : 
      71                 :            :         /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix.
      72                 :            :         *
      73                 :            :         * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
      74                 :            :         *    computed. By default, the upper triangular part is used, but can be changed
      75                 :            :         *    through the template parameter.
      76                 :            :         * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem.
      77                 :            :         * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
      78                 :            :         *    Must be less than the size of the input matrix, or an error is returned.
      79                 :            :         * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
      80                 :            :         *    respective meanings to find the largest magnitude , smallest magnitude,
      81                 :            :         *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
      82                 :            :         *    value can contain floating point value in string form, in which case the
      83                 :            :         *    eigenvalues closest to this value will be found.
      84                 :            :         * \param[in] parameters Can be ComputeEigenvectors (default) or EigenvaluesOnly.
      85                 :            :         * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
      86                 :            :         *    means machine precision.
      87                 :            :         *
      88                 :            :         * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar)
      89                 :            :         * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if
      90                 :            :         * \p parameters equals ComputeEigenvectors.
      91                 :            :         *
      92                 :            :         */
      93                 :          4 :         ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, const RMatrixType& B,
      94                 :            :                                                 Index nbrEigenvalues, std::string eigs_sigma="LM",
      95                 :            :                                                 int parameters=ComputeEigenvectors, RealScalar tol=0.0) : 
      96                 :            :                 m_eivec(), m_eivalues(), m_info(), m_isInitialized(false), m_eigenvectorsOk(false),
      97 [ +  - ][ +  - ]:          4 :                 m_nbrConverged(0), m_nbrIterations(0)
      98                 :            :         {
      99 [ +  - ][ +  - ]:          4 :                 compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     100                 :          4 :         }
     101                 :            : 
     102                 :            :         /** \brief Constructor; computes eigenvalues of given matrix.
     103                 :            :         *
     104                 :            :         * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will
     105                 :            :         *    computed. By default, the upper triangular part is used, but can be changed
     106                 :            :         *    through the template parameter.
     107                 :            :         * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
     108                 :            :         *    Must be less than the size of the input matrix, or an error is returned.
     109                 :            :         * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
     110                 :            :         *    respective meanings to find the largest magnitude , smallest magnitude,
     111                 :            :         *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
     112                 :            :         *    value can contain floating point value in string form, in which case the
     113                 :            :         *    eigenvalues closest to this value will be found.
     114                 :            :         * \param[in] parameters Can be ComputeEigenvectors (default) or EigenvaluesOnly.
     115                 :            :         * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
     116                 :            :         *    means machine precision.
     117                 :            :         *
     118                 :            :         * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar)
     119                 :            :         * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if
     120                 :            :         * \p parameters equals ComputeEigenvectors.
     121                 :            :         *
     122                 :            :         */
     123                 :         23 :         ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, Index nbrEigenvalues, std::string eigs_sigma="LM",
     124                 :            :                                                 int parameters=ComputeEigenvectors, RealScalar tol=0.0) : 
     125                 :            :                 m_eivec(), m_eivalues(), m_info(), m_isInitialized(false), m_eigenvectorsOk(false),
     126 [ +  - ][ +  - ]:         23 :                 m_nbrConverged(0), m_nbrIterations(0)
         [ +  - ][ +  - ]
     127                 :            :         {
     128 [ +  - ][ +  - ]:         23 :                 compute(A, nbrEigenvalues, eigs_sigma, parameters, tol);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     129                 :         23 :         }
     130                 :            : 
     131                 :            :         /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library.
     132                 :            :         *
     133                 :            :         * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
     134                 :            :         * \param[in]  B  Selfadjoint matrix for generalized eigenvalues.
     135                 :            :         * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
     136                 :            :         *    Must be less than the size of the input matrix, or an error is returned.
     137                 :            :         * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
     138                 :            :         *    respective meanings to find the largest magnitude , smallest magnitude,
     139                 :            :         *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
     140                 :            :         *    value can contain floating point value in string form, in which case the
     141                 :            :         *    eigenvalues closest to this value will be found.
     142                 :            :         * \param[in] parameters Can be ComputeEigenvectors (default) or EigenvaluesOnly.
     143                 :            :         * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
     144                 :            :         *    means machine precision.
     145                 :            :         *
     146                 :            :         * \returns    Reference to \c *this
     147                 :            :         *
     148                 :            :         * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK.  The eigenvalues()
     149                 :            :         * function can be used to retrieve them.  If \p parameters equals ComputeEigenvectors,
     150                 :            :         * then the eigenvectors are also computed and can be retrieved by
     151                 :            :         * calling eigenvectors().
     152                 :            :         *
     153                 :            :         */
     154                 :            :         ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, const RMatrixType& B,
     155                 :            :                                                          Index nbrEigenvalues, std::string eigs_sigma="LM",
     156                 :            :                                                          int parameters=ComputeEigenvectors, RealScalar tol=0.0);
     157                 :            : 
     158                 :            :         /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library.
     159                 :            :         *
     160                 :            :         * \param[in]  A  Selfadjoint matrix whose eigendecomposition is to be computed.
     161                 :            :         * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute.
     162                 :            :         *    Must be less than the size of the input matrix, or an error is returned.
     163                 :            :         * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with
     164                 :            :         *    respective meanings to find the largest magnitude , smallest magnitude,
     165                 :            :         *    largest algebraic, or smallest algebraic eigenvalues. Alternatively, this
     166                 :            :         *    value can contain floating point value in string form, in which case the
     167                 :            :         *    eigenvalues closest to this value will be found.
     168                 :            :         * \param[in] parameters Can be ComputeEigenvectors (default) or EigenvaluesOnly.
     169                 :            :         * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which
     170                 :            :         *    means machine precision.
     171                 :            :         *
     172                 :            :         * \returns    Reference to \c *this
     173                 :            :         *
     174                 :            :         * This function computes the eigenvalues of \p A using ARPACK.  The eigenvalues()
     175                 :            :         * function can be used to retrieve them.  If \p parameters equals ComputeEigenvectors,
     176                 :            :         * then the eigenvectors are also computed and can be retrieved by
     177                 :            :         * calling eigenvectors().
     178                 :            :         *
     179                 :            :         */
     180                 :            :         ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, Index nbrEigenvalues, 
     181                 :            :                                                          std::string eigs_sigma="LM",
     182                 :            :                                                          int parameters=ComputeEigenvectors, RealScalar tol=0.0);
     183                 :            : 
     184                 :            :         /** \brief Returns the eigenvectors of given matrix.
     185                 :            :         *
     186                 :            :         * \returns  A const reference to the matrix whose columns are the eigenvectors.
     187                 :            :         *
     188                 :            :         * \pre The eigenvectors have been computed before.
     189                 :            :         *
     190                 :            :         * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
     191                 :            :         * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
     192                 :            :         * eigenvectors are normalized to have (Euclidean) norm equal to one. If
     193                 :            :         * this object was used to solve the eigenproblem for the selfadjoint
     194                 :            :         * matrix \f$ A \f$, then the matrix returned by this function is the
     195                 :            :         * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$.
     196                 :            :         * For the generalized eigenproblem, the matrix returned is the solution of \f$A V = D B V \f$
     197                 :            :         *
     198                 :            :         * \sa eigenvalues()
     199                 :            :         */
     200                 :         27 :         const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const
     201                 :            :         {
     202 [ -  + ][ -  + ]:         27 :                 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     203 [ -  + ][ -  + ]:         27 :                 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     204                 :         27 :                 return m_eivec;
     205                 :            :         }
     206                 :            : 
     207                 :            :         /** \brief Returns the eigenvalues of given matrix.
     208                 :            :         *
     209                 :            :         * \returns A const reference to the column vector containing the eigenvalues.
     210                 :            :         *
     211                 :            :         * \pre The eigenvalues have been computed before.
     212                 :            :         *
     213                 :            :         * The eigenvalues are repeated according to their algebraic multiplicity,
     214                 :            :         * so there are as many eigenvalues as rows in the matrix. The eigenvalues
     215                 :            :         * are sorted in increasing order.
     216                 :            :         *
     217                 :            :         * \sa eigenvectors(), MatrixBase::eigenvalues()
     218                 :            :         */
     219                 :         27 :         const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
     220                 :            :         {
     221 [ -  + ][ -  + ]:         27 :                 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     222                 :         27 :                 return m_eivalues;
     223                 :            :         }
     224                 :            : 
     225                 :            :         /** \brief Computes the positive-definite square root of the matrix.
     226                 :            :         *
     227                 :            :         * \returns the positive-definite square root of the matrix
     228                 :            :         *
     229                 :            :         * \pre The eigenvalues and eigenvectors of a positive-definite matrix
     230                 :            :         * have been computed before.
     231                 :            :         *
     232                 :            :         * The square root of a positive-definite matrix \f$ A \f$ is the
     233                 :            :         * positive-definite matrix whose square equals \f$ A \f$. This function
     234                 :            :         * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
     235                 :            :         * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
     236                 :            :         *
     237                 :            :         * \sa operatorInverseSqrt(),
     238                 :            :         */
     239                 :            :         Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const
     240                 :            :         {
     241                 :            :                 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
     242                 :            :                 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
     243                 :            :                 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     244                 :            :         }
     245                 :            : 
     246                 :            :         /** \brief Computes the inverse square root of the matrix.
     247                 :            :         *
     248                 :            :         * \returns the inverse positive-definite square root of the matrix
     249                 :            :         *
     250                 :            :         * \pre The eigenvalues and eigenvectors of a positive-definite matrix
     251                 :            :         * have been computed before.
     252                 :            :         *
     253                 :            :         * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
     254                 :            :         * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
     255                 :            :         * cheaper than first computing the square root with operatorSqrt() and
     256                 :            :         * then its inverse with MatrixBase::inverse().
     257                 :            :         *
     258                 :            :         * \sa operatorSqrt(), MatrixBase::inverse(),
     259                 :            :         */
     260                 :            :         Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const
     261                 :            :         {
     262                 :            :                 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
     263                 :            :                 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
     264                 :            :                 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     265                 :            :         }
     266                 :            : 
     267                 :            :         /** \brief Reports whether previous computation was successful.
     268                 :            :         *
     269                 :            :         * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
     270                 :            :         */
     271                 :         27 :         ComputationInfo info() const
     272                 :            :         {
     273 [ -  + ][ -  + ]:         27 :                 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     274                 :         27 :                 return m_info;
     275                 :            :         }
     276                 :            : 
     277                 :            :         size_t getNbrConvergedEigenValues() const
     278                 :            :         {
     279                 :            :                 return m_nbrConverged;
     280                 :            :         }
     281                 :            : 
     282                 :         27 :         size_t getNbrIterations() const
     283                 :            :         {
     284                 :         27 :                 return m_nbrIterations; 
     285                 :            :         }
     286                 :            : 
     287                 :            : protected:
     288                 :            : 
     289                 :            :         Matrix<Scalar, Dynamic, Dynamic> m_eivec;
     290                 :            :         Matrix<Scalar, Dynamic, 1> m_eivalues;
     291                 :            :         ComputationInfo m_info;
     292                 :            :         bool m_isInitialized;
     293                 :            :         bool m_eigenvectorsOk;
     294                 :            : 
     295                 :            :         size_t m_nbrConverged;
     296                 :            :         size_t m_nbrIterations;
     297                 :            : };
     298                 :            : 
     299                 :            : template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
     300                 :            : ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>&
     301                 :         23 :     ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>
     302                 :            : ::compute(const LMatrixType& A, Index nbrEigenvalues,
     303                 :            :           std::string eigs_sigma, int parameters, RealScalar tol)
     304                 :            : {
     305 [ +  - ][ +  - ]:         23 :         RMatrixType B(0,0);
         [ +  - ][ +  - ]
     306 [ +  - ][ +  - ]:         23 :         return compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     307                 :            : }
     308                 :            : 
     309                 :            : template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
     310                 :            : ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>&
     311                 :         27 :     ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>
     312                 :            : ::compute(const LMatrixType& A, const RMatrixType& B, Index nbrEigenvalues,
     313                 :            :           std::string eigs_sigma, int parameters, RealScalar tol)
     314                 :            : {
     315 [ -  + ][ -  + ]:         27 :         eigen_assert(B.cols() == B.rows());
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     316 [ -  + ][ #  # ]:         27 :         eigen_assert(B.rows() == 0 || A.rows() == B.rows() || A.cols() == B.cols());
         [ #  # ][ -  + ]
         [ #  # ][ #  # ]
         [ -  + ][ #  # ]
         [ #  # ][ -  + ]
         [ #  # ][ #  # ]
         [ +  - ][ -  + ]
         [ #  # ][ +  - ]
         [ -  + ][ #  # ]
     317 [ +  - ][ -  + ]:         27 :         eigen_assert((parameters &~ (Eigen::EigVecMask | Eigen::GenEigMask)) == 0
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
     318                 :            :                      && (parameters & Eigen::EigVecMask) != Eigen::EigVecMask
     319                 :            :                      && "invalid option parameter");
     320                 :            : 
     321 [ -  + ][ #  # ]:         27 :         bool isBempty = (B.rows() == 0) || (B.cols() == 0);
         [ -  + ][ #  # ]
         [ -  + ][ #  # ]
         [ -  + ][ #  # ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
     322                 :            : 
     323                 :            :         // For clarity, all parameters match their ARPACK name
     324                 :            :         // Always 0 on the first call
     325                 :         27 :         int ido = 0;
     326                 :         27 :         int n = (int)A.rows();
     327                 :            :         // User parameters: "LA", "SA", "SM", "LM", "BE"
     328                 :         27 :         char whch[3] = "LM";
     329                 :            : 
     330                 :            :         // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
     331                 :         27 :         RealScalar sigma = 0.0;
     332                 :            : 
     333 [ +  - ][ +  - ]:         27 :         if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
                 [ +  - ]
           [ +  -  +  - ]
         [ +  - ][ +  - ]
           [ +  -  +  - ]
         [ +  - ][ +  - ]
           [ +  -  +  - ]
         [ +  - ][ +  - ]
           [ +  -  +  - ]
         [ +  - ][ +  - ]
           [ +  -  +  - ]
         [ +  - ][ +  - ]
                 [ +  - ]
     334                 :            :         {
     335                 :         27 :                 eigs_sigma[0] = toupper(eigs_sigma[0]);
     336                 :         27 :                 eigs_sigma[1] = toupper(eigs_sigma[1]);
     337                 :            :                 // In the following special case we're going to invert the problem, since solving
     338                 :            :                 // for larger magnitude is much much faster
     339                 :            :                 // i.e., if 'SM' is specified, we're going to really use 'LM', the default
     340 [ +  - ][ +  + ]:         27 :                 if (eigs_sigma.substr(0,2) != "SM")
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
     341                 :            :                 {
     342                 :         19 :                         whch[0] = eigs_sigma[0];
     343                 :         19 :                         whch[1] = eigs_sigma[1];
     344                 :            :                 }
     345                 :            :         }
     346                 :            :         else
     347                 :            :         {
     348                 :          0 :                 eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
     349                 :            :                 // If it's not scalar values, then the user may be explicitly
     350                 :            :                 // specifying the sigma value to cluster the evs around
     351                 :            :                 sigma = atof(eigs_sigma.c_str());
     352                 :            :                 // If atof fails, it returns 0.0, which is a fine default
     353                 :            :         }
     354                 :            : 
     355                 :            :         // "I" means normal eigenvalue problem, "G" means generalized
     356                 :            :         //
     357                 :         27 :         char bmat[2] = "I";
     358 [ +  - ][ +  - ]:         27 :         if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
         [ +  + ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ -  + ]
         [ +  - ][ +  + ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ -  + ][ +  - ]
         [ -  + ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
         [ #  # ][ +  - ]
         [ +  - ][ -  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ #  # ]
         [ +  - ][ +  - ]
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ #  # ][ +  - ]
         [ +  - ][ -  + ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ #  # ]
     359                 :          8 :                 bmat[0] = 'G';
     360                 :            : 
     361                 :            :         // Now we determine the mode to use
     362                 :            :         //
     363 [ +  + ][ -  + ]:         27 :         int mode = (bmat[0] == 'G') + 1;
         [ -  + ][ +  - ]
         [ +  - ][ +  - ]
     364 [ +  - ][ +  - ]:         27 :         if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
         [ +  + ][ +  - ]
         [ +  - ][ +  - ]
         [ -  + ][ +  - ]
         [ +  + ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ -  + ][ +  - ]
         [ -  + ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ -  + ][ +  - ]
         [ -  + ][ #  # ]
         [ +  - ][ +  - ]
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ #  # ]
         [ +  - ][ +  - ]
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ #  # ]
         [ +  - ][ +  - ]
         [ -  + ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ #  # ]
     365                 :            :         {
     366                 :            :                 // We're going to use shift-and-invert mode, and basically find
     367                 :            :                 // the largest eigenvalues of the inverse operator
     368                 :          8 :                 mode = 3;
     369                 :            :                 //std::cout << "Mode = 3\n";
     370                 :            :         }
     371                 :            : 
     372                 :            :         // The user-specified number of eigenvalues/vectors to compute
     373                 :         27 :         int nev = (int)nbrEigenvalues;
     374                 :            :         // Allocate space for ARPACK to store the residual
     375                 :         27 :         Scalar *resid = new Scalar[n];
     376                 :            :         // Number of Lanczos vectors, must satisfy nev < ncv <= n
     377                 :            :         // Note that this indicates that nev != n, and we cannot compute
     378                 :            :         // all eigenvalues of a mtrix
     379                 :         27 :         int ncv = std::min(std::max(2*nev, 20), n);
     380                 :            : 
     381                 :            :         // The working n x ncv matrix, also store the final eigenvectors (if computed)
     382                 :         27 :         Scalar *v = new Scalar[n*ncv];
     383                 :         27 :         int ldv = n;
     384                 :            : 
     385                 :            :         // Working space
     386                 :         27 :         Scalar *workd = new Scalar[3*n];
     387                 :         27 :         int lworkl = ncv*ncv+8*ncv; // Must be at least this length
     388                 :         27 :         Scalar *workl = new Scalar[lworkl];
     389                 :            : 
     390 [ +  + ][ +  + ]:        324 :         int *iparam= new int[11];
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     391                 :         27 :         iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
     392                 :         27 :         iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1)));
     393                 :         27 :         iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
     394                 :            : 
     395                 :            :         // Used during reverse communicate to notify where arrays start
     396 [ +  + ][ +  + ]:        324 :         int *ipntr = new int[11]; 
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     397                 :            : 
     398                 :            :         // Error codes are returned in here, initial value of 0 indicates a random initial
     399                 :            :         // residual vector is used, any other values means resid contains the initial residual
     400                 :            :         // vector, possibly from a previous run
     401                 :         27 :         int cinfo = 0;
     402                 :            : 
     403                 :         27 :         Scalar scale = 1.0;
     404                 :            :         /*
     405                 :            :         if (!isBempty)
     406                 :            :         {
     407                 :            :         Scalar scale = B.norm() / std::sqrt(n);
     408                 :            :         scale = std::pow(2, std::floor(std::log(scale+1)));
     409                 :            :         ////M /= scale;
     410                 :            :         for (size_t i=0; i<(size_t)B.outerSize(); i++)
     411                 :            :             for (typename RMatrixType::InnerIterator it(B, i); it; ++it)
     412                 :            :                 it.valueRef() /= scale;
     413                 :            :         }
     414                 :            :         */
     415                 :            : 
     416                 :         27 :         MatrixOperation op(A);
     417                 :            :         //              (mode==1 || mode==2) ? B :
     418                 :            :         //                   (mode==3) ? A : LMatrixType());
     419                 :            : 
     420 [ +  + ][ +  + ]:        761 :         do
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     421                 :            :         {
     422                 :            :                 //std::cout << "Entering main loop\n";
     423 [ +  - ][ +  - ]:        761 :                 arpack_internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, 
         [ +  - ][ +  - ]
     424                 :            :                                                                            &ncv, v, &ldv, iparam, ipntr, workd, workl,
     425                 :            :                                                                            &lworkl, &cinfo);
     426         [ +  - ]:        761 :                 if (ido == -1 || ido == 1)
           [ +  +  +  - ]
           [ +  +  +  - ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
           [ +  +  +  - ]
                 [ +  + ]
     427                 :            :                 {
     428                 :        445 :                         Scalar *in  = workd + ipntr[0] - 1;
     429                 :        445 :                         Scalar *out = workd + ipntr[1] - 1;
     430                 :            : 
     431 [ +  + ][ +  - ]:        445 :                         if (ido == 1 && mode != 2)
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  + ][ +  - ]
         [ +  + ][ +  - ]
         [ +  + ][ +  - ]
     432                 :            :                         {
     433                 :        437 :                                 Scalar *out2 = workd + ipntr[2] - 1;
     434 [ -  + ][ #  # ]:        437 :                                 if (isBempty || mode == 1)
         [ -  + ][ #  # ]
         [ -  + ][ #  # ]
         [ -  + ][ #  # ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
     435 [ +  - ][ +  - ]:        408 :                                         Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
     436                 :            :                                 else
     437 [ #  # ][ #  # ]:         29 :                                         Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
            [ # ][ #  # ]
                 [ #  # ]
     438                 :        437 :                                 in = workd + ipntr[2] - 1;
     439                 :            :                         }
     440                 :            : 
     441 [ +  + ][ +  - ]:        445 :                         if (mode == 1)
         [ +  - ][ -  + ]
         [ -  + ][ -  + ]
     442                 :            :                         {
     443 [ +  - ][ +  - ]:        345 :                                 if (isBempty)
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
     444                 :            :                                 {
     445                 :            :                                         // OP = A
     446 [ +  - ][ +  - ]:        345 :                                         Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
               [ + ][ + ]
         [ +  - ][ +  - ]
     447                 :            :                                 }
     448                 :            :                                 else
     449                 :            :                                 {
     450                 :            :                                         // TODO OP = L^{-1}AL^{-T}
     451                 :            :                                         //internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
     452                 :            :                                 }
     453                 :            :                         }
     454 [ -  + ][ #  # ]:        100 :                         else if (mode == 2)
         [ #  # ][ -  + ]
         [ -  + ][ -  + ]
     455                 :            :                         {
     456 [ #  # ][ #  # ]:          0 :                                 if (ido == 1)
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     457 [ #  # ][ #  # ]:          0 :                                         Matrix<Scalar, Dynamic, 1>::Map(in, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
               [ # ][ # ]
         [ #  # ][ #  # ]
     458                 :            :                                 // OP = B^{-1} A
     459 [ #  # ][ #  # ]:          0 :                                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
               [ # ][ # ]
         [ #  # ][ #  # ]
     460                 :            :                         }
     461 [ +  - ][ #  # ]:        100 :                         else if (mode == 3)
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
     462                 :            :                         {
     463                 :            :                                 // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
     464                 :            :                                 // The B * in is already computed and stored at in if ido == 1
     465 [ +  + ][ +  - ]:        100 :                                 if (ido == 1 || isBempty)
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  + ][ +  - ]
         [ +  + ][ -  + ]
         [ +  + ][ -  + ]
     466 [ #  # ][ #  # ]:         96 :                                         Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
     467                 :            :                                 else
     468 [ #  # ][ #  # ]:        445 :                                         Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
            [ #  # ][ # ]
          [ # ][ # ][ # ]
         [ #  # ][ #  # ]
     469                 :            :                         }
     470                 :            :                 }
     471 [ +  + ][ -  + ]:        316 :                 else if (ido == 2)
         [ -  + ][ +  + ]
         [ +  + ][ +  + ]
     472                 :            :                 {
     473                 :        289 :                         Scalar *in  = workd + ipntr[0] - 1;
     474                 :        289 :                         Scalar *out = workd + ipntr[1] - 1;
     475                 :            : 
     476 [ -  + ][ #  # ]:        289 :                         if (isBempty || mode == 1)
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ -  + ][ #  # ]
         [ +  - ][ -  + ]
         [ +  - ][ -  + ]
     477 [ +  - ][ +  - ]:        194 :                                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
     478                 :            :                         else
     479 [ #  # ][ #  # ]:        289 :                                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     480                 :            :                 }
     481                 :            :                 //Scalar *out = workd + ipntr[1] - 1;
     482                 :            :                 //std::cout << Matrix<Scalar, Dynamic, 1>::Map(out, n).transpose() << std::endl;
     483                 :            :         } while (ido != 99);
     484                 :            :                 
     485                 :            :         //std::cout << "Finsihed\n";
     486                 :            : 
     487 [ -  + ][ -  + ]:         27 :         if (cinfo == 1)
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     488                 :            :         {
     489                 :            :                 //std::cout << "FAILED WITH NO CONV\n";
     490                 :          0 :                 m_info = NoConvergence;
     491                 :            :         }
     492 [ -  + ][ -  + ]:         27 :         else if (cinfo == 3)
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     493                 :            :         {
     494                 :            :                 //std::cout << "FAILED WITH NUMERICAL ISSUE\n";
     495                 :          0 :                 m_info = NumericalIssue;
     496                 :            :         }
     497 [ -  + ][ -  + ]:         27 :         else if (cinfo < 0)
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     498                 :            :         {
     499                 :            :                 //std::cout << "FAILED WITH INVALID INPUT " << info << "\n";
     500                 :            : #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
     501                 :            :                 m_info = NoConvergence;
     502                 :            : #else
     503                 :          0 :                 m_info = InvalidInput;
     504                 :            : #endif
     505                 :            :         }
     506 [ -  + ][ -  + ]:         27 :         else if (cinfo != 0)
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     507                 :          0 :                 eigen_assert(false && "Unknown ARPACK return value!");
     508                 :            :         else
     509                 :            :         {
     510                 :            :                 // Do we compute eigenvectors or not?
     511                 :         27 :                 int rvec = (parameters & ComputeEigenvectors) == ComputeEigenvectors;
     512                 :            : 
     513                 :            :                 // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
     514                 :         27 :                 char howmny[2] = "A"; 
     515                 :            : 
     516                 :            :                 // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
     517 [ +  + ][ +  + ]:        464 :                 int *select = new int[ncv];
         [ +  + ][ +  - ]
         [ +  + ][ +  - ]
         [ +  + ][ +  - ]
         [ +  + ][ +  - ]
         [ +  + ][ +  + ]
     518                 :            : 
     519                 :            :                 // Final eigenvalues
     520 [ +  - ][ +  - ]:         27 :                 m_eivalues.resize(nev, 1);
         [ +  - ][ +  - ]
     521                 :            : 
     522 [ +  - ][ +  - ]:         27 :                 arpack_internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
         [ +  - ][ +  - ]
     523                 :            :                                                                            &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
     524                 :            :                                                                            v, &ldv, iparam, ipntr, workd, workl, &lworkl, &cinfo);
     525                 :            : 
     526   [ -  +  -  +  :         27 :                 if (cinfo == -14)
           -  + ][ -  + ]
         [ -  + ][ -  + ]
           [ -  +  -  + ]
     527                 :          0 :                         m_info = NoConvergence;
     528 [ -  + ][ -  + ]:         27 :                 else if (cinfo != 0)
         [ -  + ][ -  + ]
         [ -  + ][ -  + ]
     529                 :            : #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
     530                 :            :                         m_info = NoConvergence;
     531                 :            : #else
     532                 :          0 :                         m_info = InvalidInput;
     533                 :            : #endif
     534                 :            :                 else
     535                 :            :                 {
     536 [ +  - ][ +  - ]:         27 :                         if (rvec)
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     537                 :            :                         {
     538 [ +  - ][ +  - ]:         27 :                                 m_eivec.resize(A.rows(), nev);
         [ +  - ][ +  - ]
                 [ +  - ]
     539 [ +  + ][ +  + ]:         83 :                                 for (int i=0; i<nev; i++)
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     540 [ +  + ][ +  + ]:       1625 :                                         for (int j=0; j<n; j++)
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     541 [ +  - ][ +  - ]:       1569 :                                                 m_eivec(j,i) = v[i*n+j] / scale;
         [ +  - ][ +  - ]
     542                 :            :   
     543                 :            :                                 // TODO if (mode == 1 && !isBempty && BisSPD)
     544                 :            :                                 //  internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
     545                 :            : 
     546                 :         27 :                                 m_eigenvectorsOk = true;
     547                 :            :                         }
     548                 :            : 
     549                 :         27 :                         m_nbrIterations = iparam[2];
     550                 :         27 :                         m_nbrConverged  = iparam[4];
     551                 :         27 :                         m_info = Success;
     552                 :            :                 }
     553 [ +  - ][ +  - ]:         27 :                 delete[] select;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     554                 :            :         }
     555                 :            : 
     556 [ +  - ][ +  - ]:         27 :         delete[] v;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     557 [ +  - ][ +  - ]:         27 :         delete[] iparam;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     558 [ +  - ][ +  - ]:         27 :         delete[] ipntr;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     559 [ +  - ][ +  - ]:         27 :         delete[] workd;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     560 [ +  - ][ +  - ]:         27 :         delete[] workl;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     561 [ +  - ][ +  - ]:         27 :         delete[] resid;
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
     562                 :         27 :         m_isInitialized = true;
     563                 :         27 :         return *this;
     564                 :            : }
     565                 :            : 
     566                 :            : // Single precision
     567                 :            : //
     568                 :            : extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which,
     569                 :            :                         int *nev, float *tol, float *resid, int *ncv,
     570                 :            :                         float *v, int *ldv, int *iparam, int *ipntr,
     571                 :            :                         float *workd, float *workl, int *lworkl,
     572                 :            :                         int *info);
     573                 :            : 
     574                 :            : extern "C" void sseupd_(int *rvec, char *All, int *select, float *d,
     575                 :            :                         float *z, int *ldz, float *sigma, 
     576                 :            :                         char *bmat, int *n, char *which, int *nev,
     577                 :            :                         float *tol, float *resid, int *ncv, float *v,
     578                 :            :                         int *ldv, int *iparam, int *ipntr, float *workd,
     579                 :            :                         float *workl, int *lworkl, int *ierr);
     580                 :            : 
     581                 :            : // Double precision
     582                 :            : //
     583                 :            : extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
     584                 :            :                         int *nev, double *tol, double *resid, int *ncv,
     585                 :            :                         double *v, int *ldv, int *iparam, int *ipntr,
     586                 :            :                         double *workd, double *workl, int *lworkl,
     587                 :            :                         int *info);
     588                 :            : 
     589                 :            : extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
     590                 :            :                         double *z, int *ldz, double *sigma, 
     591                 :            :                         char *bmat, int *n, char *which, int *nev,
     592                 :            :                         double *tol, double *resid, int *ncv, double *v,
     593                 :            :                         int *ldv, int *iparam, int *ipntr, double *workd,
     594                 :            :                         double *workl, int *lworkl, int *ierr);
     595                 :            : 
     596                 :            : namespace arpack_internal {
     597                 :            : 
     598                 :            : template<typename Scalar, typename RealScalar> struct arpack_wrapper
     599                 :            : {
     600                 :            :   static inline void saupd(int *ido, char *bmat, int *n, char *which,
     601                 :            :       int *nev, RealScalar *tol, Scalar *resid, int *ncv,
     602                 :            :       Scalar *v, int *ldv, int *iparam, int *ipntr,
     603                 :            :       Scalar *workd, Scalar *workl, int *lworkl, int *info);
     604                 :            : 
     605                 :            :   static inline void seupd(int *rvec, char *All, int *select, Scalar *d,
     606                 :            :       Scalar *z, int *ldz, RealScalar *sigma,
     607                 :            :       char *bmat, int *n, char *which, int *nev,
     608                 :            :       RealScalar *tol, Scalar *resid, int *ncv, Scalar *v,
     609                 :            :       int *ldv, int *iparam, int *ipntr, Scalar *workd,
     610                 :            :       Scalar *workl, int *lworkl, int *ierr);
     611                 :            : };
     612                 :            : 
     613                 :            : template <> struct arpack_wrapper<float, float>
     614                 :            : {
     615                 :            :   static inline void saupd(int *ido, char *bmat, int *n, char *which,
     616                 :            :       int *nev, float *tol, float *resid, int *ncv,
     617                 :            :       float *v, int *ldv, int *iparam, int *ipntr,
     618                 :            :       float *workd, float *workl, int *lworkl, int *info)
     619                 :            :   {
     620                 :            :     ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
     621                 :            :   }
     622                 :            : 
     623                 :            :   static inline void seupd(int *rvec, char *All, int *select, float *d,
     624                 :            :       float* z, int* ldz, float *sigma,
     625                 :            :       char *bmat, int *n, char *which, int *nev,
     626                 :            :       float *tol, float *resid, int *ncv, float *v,
     627                 :            :       int *ldv, int *iparam, int *ipntr, float *workd,
     628                 :            :       float *workl, int *lworkl, int *ierr)
     629                 :            :   {
     630                 :            :     sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
     631                 :            :         workd, workl, lworkl, ierr);
     632                 :            :   }
     633                 :            : };
     634                 :            : 
     635                 :            : template <> struct arpack_wrapper<double, double>
     636                 :            : {
     637                 :        761 :   static inline void saupd(int *ido, char *bmat, int *n, char *which,
     638                 :            :       int *nev, double *tol, double *resid, int *ncv,
     639                 :            :       double *v, int *ldv, int *iparam, int *ipntr,
     640                 :            :       double *workd, double *workl, int *lworkl, int *info)
     641                 :            :   {
     642                 :        761 :     dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
     643                 :        761 :   }
     644                 :            : 
     645                 :         27 :   static inline void seupd(int *rvec, char *All, int *select, double *d,
     646                 :            :       double* z, int* ldz, double *sigma,
     647                 :            :       char *bmat, int *n, char *which, int *nev,
     648                 :            :       double *tol, double *resid, int *ncv, double *v,
     649                 :            :       int *ldv, int *iparam, int *ipntr, double *workd,
     650                 :            :       double *workl, int *lworkl, int *ierr)
     651                 :            :   {
     652                 :            :     dseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
     653                 :         27 :         workd, workl, lworkl, ierr);
     654                 :         27 :   }
     655                 :            : };
     656                 :            : 
     657                 :            : }
     658                 :            : 
     659                 :            : #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
     660                 :            : 

Generated by: LCOV version 1.9