Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
arpack_wrapper.hpp
Go to the documentation of this file.
1 /* This software is distributed under BSD 3-clause license (see LICENSE file).
2  *
3  * This code is based on Eigen3 ARPACK wrapper
4  * Copyright (c) 2012 David Harmon
5  *
6  * The code is relicensed under BSD 3-clause with the permission of its
7  * original author. Some changes were made afterwards by Sergey Lisitsyn.
8  */
9 
10 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
11 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
12 
13 /* Tapkee includes */
14 #include <tapkee/defines.hpp>
15 /* End of Tapkee includes */
16 
17 using Eigen::Matrix;
18 using Eigen::Dynamic;
19 using Eigen::ComputationInfo;
20 using Eigen::ComputeEigenvectors;
21 using Eigen::NoConvergence;
22 using Eigen::NumericalIssue;
23 #ifndef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
24 using Eigen::InvalidInput;
25 #endif
26 using Eigen::Success;
27 using Eigen::NumTraits;
28 
30 namespace arpack_internal
31 {
32  template<typename Scalar, typename RealScalar> struct arpack_wrapper;
33 }
34 
35 template<class LMatrixType, class RMatrixType, class MatrixOperation, bool BisSPD=false>
37 {
38 public:
39  //typedef typename MatrixSolver::MatrixType MatrixType;
40 
42  typedef typename LMatrixType::Scalar Scalar;
43  typedef typename LMatrixType::Index Index;
44 
51  typedef typename NumTraits<Scalar>::Real RealScalar;
52 
58  typedef typename Eigen::internal::plain_col_type<LMatrixType, RealScalar>::type RealVectorType;
59 
68  {
69  }
70 
93  ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, const RMatrixType& B,
94  Index nbrEigenvalues, std::string eigs_sigma="LM",
95  int parameters=ComputeEigenvectors, RealScalar tol=0.0) :
98  {
99  compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
100  }
101 
123  ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, Index nbrEigenvalues, std::string eigs_sigma="LM",
124  int parameters=ComputeEigenvectors, RealScalar tol=0.0) :
125  m_eivec(), m_eivalues(), m_info(), m_isInitialized(false), m_eigenvectorsOk(false),
127  {
128  compute(A, nbrEigenvalues, eigs_sigma, parameters, tol);
129  }
130 
154  ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, const RMatrixType& B,
155  Index nbrEigenvalues, std::string eigs_sigma="LM",
156  int parameters=ComputeEigenvectors, RealScalar tol=0.0);
157 
180  ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, Index nbrEigenvalues,
181  std::string eigs_sigma="LM",
182  int parameters=ComputeEigenvectors, RealScalar tol=0.0);
183 
200  const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const
201  {
202  eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
203  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
204  return m_eivec;
205  }
206 
219  const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
220  {
221  eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
222  return m_eivalues;
223  }
224 
239  Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const
240  {
241  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
242  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
243  return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
244  }
245 
260  Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const
261  {
262  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
263  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
264  return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
265  }
266 
271  ComputationInfo info() const
272  {
273  eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
274  return m_info;
275  }
276 
278  {
279  return m_nbrConverged;
280  }
281 
282  size_t getNbrIterations() const
283  {
284  return m_nbrIterations;
285  }
286 
287 protected:
288 
289  Matrix<Scalar, Dynamic, Dynamic> m_eivec;
290  Matrix<Scalar, Dynamic, 1> m_eivalues;
291  ComputationInfo m_info;
294 
297 };
298 
299 template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
302 ::compute(const LMatrixType& A, Index nbrEigenvalues,
303  std::string eigs_sigma, int parameters, RealScalar tol)
304 {
305  RMatrixType B(0,0);
306  return compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
307 }
308 
309 template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
312 ::compute(const LMatrixType& A, const RMatrixType& B, Index nbrEigenvalues,
313  std::string eigs_sigma, int parameters, RealScalar tol)
314 {
315  eigen_assert(B.cols() == B.rows());
316  eigen_assert(B.rows() == 0 || A.rows() == B.rows() || A.cols() == B.cols());
317  eigen_assert((parameters &~ (Eigen::EigVecMask | Eigen::GenEigMask)) == 0
318  && (parameters & Eigen::EigVecMask) != Eigen::EigVecMask
319  && "invalid option parameter");
320 
321  bool isBempty = (B.rows() == 0) || (B.cols() == 0);
322 
323  // For clarity, all parameters match their ARPACK name
324  // Always 0 on the first call
325  int ido = 0;
326  int n = (int)A.rows();
327  // User parameters: "LA", "SA", "SM", "LM", "BE"
328  char whch[3] = "LM";
329 
330  // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
331  RealScalar sigma = 0.0;
332 
333  if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
334  {
335  eigs_sigma[0] = toupper(eigs_sigma[0]);
336  eigs_sigma[1] = toupper(eigs_sigma[1]);
337  // In the following special case we're going to invert the problem, since solving
338  // for larger magnitude is much much faster
339  // i.e., if 'SM' is specified, we're going to really use 'LM', the default
340  if (eigs_sigma.substr(0,2) != "SM")
341  {
342  whch[0] = eigs_sigma[0];
343  whch[1] = eigs_sigma[1];
344  }
345  }
346  else
347  {
348  eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
349  // If it's not scalar values, then the user may be explicitly
350  // specifying the sigma value to cluster the evs around
351  sigma = atof(eigs_sigma.c_str());
352  // If atof fails, it returns 0.0, which is a fine default
353  }
354 
355  // "I" means normal eigenvalue problem, "G" means generalized
356  //
357  char bmat[2] = "I";
358  if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
359  bmat[0] = 'G';
360 
361  // Now we determine the mode to use
362  //
363  int mode = (bmat[0] == 'G') + 1;
364  if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
365  {
366  // We're going to use shift-and-invert mode, and basically find
367  // the largest eigenvalues of the inverse operator
368  mode = 3;
369  //std::cout << "Mode = 3\n";
370  }
371 
372  // The user-specified number of eigenvalues/vectors to compute
373  int nev = (int)nbrEigenvalues;
374  // Allocate space for ARPACK to store the residual
375  Scalar *resid = new Scalar[n];
376  // Number of Lanczos vectors, must satisfy nev < ncv <= n
377  // Note that this indicates that nev != n, and we cannot compute
378  // all eigenvalues of a mtrix
379  int ncv = std::min(std::max(2*nev, 20), n);
380 
381  // The working n x ncv matrix, also store the final eigenvectors (if computed)
382  Scalar *v = new Scalar[n*ncv];
383  int ldv = n;
384 
385  // Working space
386  Scalar *workd = new Scalar[3*n];
387  int lworkl = ncv*ncv+8*ncv; // Must be at least this length
388  Scalar *workl = new Scalar[lworkl];
389 
390  int *iparam= new int[11];
391  iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
392  iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1)));
393  iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
394 
395  // Used during reverse communicate to notify where arrays start
396  int *ipntr = new int[11];
397 
398  // Error codes are returned in here, initial value of 0 indicates a random initial
399  // residual vector is used, any other values means resid contains the initial residual
400  // vector, possibly from a previous run
401  int cinfo = 0;
402 
403  Scalar scale = 1.0;
404  /*
405  if (!isBempty)
406  {
407  Scalar scale = B.norm() / std::sqrt(n);
408  scale = std::pow(2, std::floor(std::log(scale+1)));
410  for (size_t i=0; i<(size_t)B.outerSize(); i++)
411  for (typename RMatrixType::InnerIterator it(B, i); it; ++it)
412  it.valueRef() /= scale;
413  }
414  */
415 
416  MatrixOperation op(A);
417  // (mode==1 || mode==2) ? B :
418  // (mode==3) ? A : LMatrixType());
419 
420  do
421  {
422  //std::cout << "Entering main loop\n";
423  arpack_internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid,
424  &ncv, v, &ldv, iparam, ipntr, workd, workl,
425  &lworkl, &cinfo);
426  if (ido == -1 || ido == 1)
427  {
428  Scalar *in = workd + ipntr[0] - 1;
429  Scalar *out = workd + ipntr[1] - 1;
430 
431  if (ido == 1 && mode != 2)
432  {
433  Scalar *out2 = workd + ipntr[2] - 1;
434  if (isBempty || mode == 1)
435  Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
436  else
437  Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
438  in = workd + ipntr[2] - 1;
439  }
440 
441  if (mode == 1)
442  {
443  if (isBempty)
444  {
445  // OP = A
446  Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
447  }
448  else
449  {
450  // TODO OP = L^{-1}AL^{-T}
451  //internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
452  }
453  }
454  else if (mode == 2)
455  {
456  if (ido == 1)
457  Matrix<Scalar, Dynamic, 1>::Map(in, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
458  // OP = B^{-1} A
459  Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
460  }
461  else if (mode == 3)
462  {
463  // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
464  // The B * in is already computed and stored at in if ido == 1
465  if (ido == 1 || isBempty)
466  Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
467  else
468  Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
469  }
470  }
471  else if (ido == 2)
472  {
473  Scalar *in = workd + ipntr[0] - 1;
474  Scalar *out = workd + ipntr[1] - 1;
475 
476  if (isBempty || mode == 1)
477  Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
478  else
479  Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
480  }
481  //Scalar *out = workd + ipntr[1] - 1;
482  //std::cout << Matrix<Scalar, Dynamic, 1>::Map(out, n).transpose() << std::endl;
483  } while (ido != 99);
484 
485  //std::cout << "Finsihed\n";
486 
487  if (cinfo == 1)
488  {
489  //std::cout << "FAILED WITH NO CONV\n";
490  m_info = NoConvergence;
491  }
492  else if (cinfo == 3)
493  {
494  //std::cout << "FAILED WITH NUMERICAL ISSUE\n";
495  m_info = NumericalIssue;
496  }
497  else if (cinfo < 0)
498  {
499  //std::cout << "FAILED WITH INVALID INPUT " << info << "\n";
500 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
501  m_info = NoConvergence;
502 #else
503  m_info = InvalidInput;
504 #endif
505  }
506  else if (cinfo != 0)
507  eigen_assert(false && "Unknown ARPACK return value!");
508  else
509  {
510  // Do we compute eigenvectors or not?
511  int rvec = (parameters & ComputeEigenvectors) == ComputeEigenvectors;
512 
513  // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
514  char howmny[2] = "A";
515 
516  // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
517  int *select = new int[ncv];
518 
519  // Final eigenvalues
520  m_eivalues.resize(nev, 1);
521 
522  arpack_internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
523  &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
524  v, &ldv, iparam, ipntr, workd, workl, &lworkl, &cinfo);
525 
526  if (cinfo == -14)
527  m_info = NoConvergence;
528  else if (cinfo != 0)
529 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
530  m_info = NoConvergence;
531 #else
532  m_info = InvalidInput;
533 #endif
534  else
535  {
536  if (rvec)
537  {
538  m_eivec.resize(A.rows(), nev);
539  for (int i=0; i<nev; i++)
540  for (int j=0; j<n; j++)
541  m_eivec(j,i) = v[i*n+j] / scale;
542 
543  // TODO if (mode == 1 && !isBempty && BisSPD)
544  // internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
545 
546  m_eigenvectorsOk = true;
547  }
548 
549  m_nbrIterations = iparam[2];
550  m_nbrConverged = iparam[4];
551  m_info = Success;
552  }
553  delete[] select;
554  }
555 
556  delete[] v;
557  delete[] iparam;
558  delete[] ipntr;
559  delete[] workd;
560  delete[] workl;
561  delete[] resid;
562  m_isInitialized = true;
563  return *this;
564 }
565 
566 // Single precision
567 //
568 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which,
569  int *nev, float *tol, float *resid, int *ncv,
570  float *v, int *ldv, int *iparam, int *ipntr,
571  float *workd, float *workl, int *lworkl,
572  int *info);
573 
574 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d,
575  float *z, int *ldz, float *sigma,
576  char *bmat, int *n, char *which, int *nev,
577  float *tol, float *resid, int *ncv, float *v,
578  int *ldv, int *iparam, int *ipntr, float *workd,
579  float *workl, int *lworkl, int *ierr);
580 
581 // Double precision
582 //
583 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
584  int *nev, double *tol, double *resid, int *ncv,
585  double *v, int *ldv, int *iparam, int *ipntr,
586  double *workd, double *workl, int *lworkl,
587  int *info);
588 
589 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
590  double *z, int *ldz, double *sigma,
591  char *bmat, int *n, char *which, int *nev,
592  double *tol, double *resid, int *ncv, double *v,
593  int *ldv, int *iparam, int *ipntr, double *workd,
594  double *workl, int *lworkl, int *ierr);
595 
596 namespace arpack_internal {
597 
598 template<typename Scalar, typename RealScalar> struct arpack_wrapper
599 {
600  static inline void saupd(int *ido, char *bmat, int *n, char *which,
601  int *nev, RealScalar *tol, Scalar *resid, int *ncv,
602  Scalar *v, int *ldv, int *iparam, int *ipntr,
603  Scalar *workd, Scalar *workl, int *lworkl, int *info);
604 
605  static inline void seupd(int *rvec, char *All, int *select, Scalar *d,
606  Scalar *z, int *ldz, RealScalar *sigma,
607  char *bmat, int *n, char *which, int *nev,
608  RealScalar *tol, Scalar *resid, int *ncv, Scalar *v,
609  int *ldv, int *iparam, int *ipntr, Scalar *workd,
610  Scalar *workl, int *lworkl, int *ierr);
611 };
612 
613 template <> struct arpack_wrapper<float, float>
614 {
615  static inline void saupd(int *ido, char *bmat, int *n, char *which,
616  int *nev, float *tol, float *resid, int *ncv,
617  float *v, int *ldv, int *iparam, int *ipntr,
618  float *workd, float *workl, int *lworkl, int *info)
619  {
620  ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
621  }
622 
623  static inline void seupd(int *rvec, char *All, int *select, float *d,
624  float* z, int* ldz, float *sigma,
625  char *bmat, int *n, char *which, int *nev,
626  float *tol, float *resid, int *ncv, float *v,
627  int *ldv, int *iparam, int *ipntr, float *workd,
628  float *workl, int *lworkl, int *ierr)
629  {
630  sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
631  workd, workl, lworkl, ierr);
632  }
633 };
634 
635 template <> struct arpack_wrapper<double, double>
636 {
637  static inline void saupd(int *ido, char *bmat, int *n, char *which,
638  int *nev, double *tol, double *resid, int *ncv,
639  double *v, int *ldv, int *iparam, int *ipntr,
640  double *workd, double *workl, int *lworkl, int *info)
641  {
642  dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
643  }
644 
645  static inline void seupd(int *rvec, char *All, int *select, double *d,
646  double* z, int* ldz, double *sigma,
647  char *bmat, int *n, char *which, int *nev,
648  double *tol, double *resid, int *ncv, double *v,
649  int *ldv, int *iparam, int *ipntr, double *workd,
650  double *workl, int *lworkl, int *ierr)
651  {
652  dseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
653  workd, workl, lworkl, ierr);
654  }
655 };
656 
657 }
658 
659 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
660 
void ssaupd_(int *ido, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *info)
void dseupd_(int *rvec, char *All, int *select, double *d, double *z, int *ldz, double *sigma, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *ierr)
ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType &A, const RMatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes generalized eigenvalues of given matrix with respect to another matrix...
void dsaupd_(int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
ArpackGeneralizedSelfAdjointEigenSolver & compute(const LMatrixType &A, const RMatrixType &B, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library...
static void seupd(int *rvec, char *All, int *select, float *d, float *z, int *ldz, float *sigma, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *ierr)
static void saupd(int *ido, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *info)
Matrix< Scalar, Dynamic, Dynamic > operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Eigen::internal::plain_col_type< LMatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
static void saupd(int *ido, char *bmat, int *n, char *which, int *nev, RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, int *ldv, int *iparam, int *ipntr, Scalar *workd, Scalar *workl, int *lworkl, int *info)
static void seupd(int *rvec, char *All, int *select, Scalar *d, Scalar *z, int *ldz, RealScalar *sigma, char *bmat, int *n, char *which, int *nev, RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, int *ldv, int *iparam, int *ipntr, Scalar *workd, Scalar *workl, int *lworkl, int *ierr)
NumTraits< Scalar >::Real RealScalar
Real scalar type for MatrixType.
ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType &A, Index nbrEigenvalues, std::string eigs_sigma="LM", int parameters=ComputeEigenvectors, RealScalar tol=0.0)
Constructor; computes eigenvalues of given matrix.
Matrix< Scalar, Dynamic, Dynamic > operatorSqrt() const
Computes the positive-definite square root of the matrix.
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
static void saupd(int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
ComputationInfo info() const
Reports whether previous computation was successful.
void sseupd_(int *rvec, char *All, int *select, float *d, float *z, int *ldz, float *sigma, char *bmat, int *n, char *which, int *nev, float *tol, float *resid, int *ncv, float *v, int *ldv, int *iparam, int *ipntr, float *workd, float *workl, int *lworkl, int *ierr)
LMatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
static void seupd(int *rvec, char *All, int *select, double *d, double *z, int *ldz, double *sigma, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *ierr)
ArpackGeneralizedSelfAdjointEigenSolver()
Default constructor.
Matrix< Scalar, Dynamic, 1 > m_eivalues
Matrix< Scalar, Dynamic, Dynamic > m_eivec