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