• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 David Harmon <dharmon@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    */
ArpackGeneralizedSelfAdjointEigenSolver()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)
m_eivec()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)
m_eivec()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    */
eigenvectors()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    */
eigenvalues()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    */
operatorSqrt()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    */
operatorInverseSqrt()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    */
info()305   ComputationInfo info() const
306   {
307     eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
308     return m_info;
309   }
310 
getNbrConvergedEigenValues()311   size_t getNbrConvergedEigenValues() const
312   { return m_nbrConverged; }
313 
getNbrIterations()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>
compute(const MatrixType & A,Index nbrEigenvalues,std::string eigs_sigma,int options,RealScalar tol)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>
compute(const MatrixType & A,const MatrixType & B,Index nbrEigenvalues,std::string eigs_sigma,int options,RealScalar tol)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 {
saupdarpack_wrapper682   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 
seupdarpack_wrapper690   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