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