Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
ArpackSelfAdjointEigenSolver.h
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 // 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 
31  typedef typename MatrixType::Scalar Scalar;
32  typedef typename MatrixType::Index Index;
33 
40  typedef typename NumTraits<Scalar>::Real RealScalar;
41 
47  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
48 
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 
86  ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B,
87  Index nbrEigenvalues, std::string eigs_sigma="LM",
88  int options=ComputeEigenvectors, RealScalar tol=0.0)
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 
121  ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A,
122  Index nbrEigenvalues, std::string eigs_sigma="LM",
123  int options=ComputeEigenvectors, RealScalar tol=0.0)
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 
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 
184  ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A,
185  Index nbrEigenvalues, std::string eigs_sigma="LM",
186  int options=ComputeEigenvectors, RealScalar tol=0.0);
187 
188 
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 
230  const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
231  {
232  eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
233  return m_eivalues;
234  }
235 
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 
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 
290  ComputationInfo info() const
291  {
292  eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
293  return m_info;
294  }
295 
296  size_t getNbrConvergedEigenValues() const
297  { return m_nbrConverged; }
298 
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>
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>
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)));
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 {
667  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 
675  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 
Namespace containing all symbols from the Eigen library.
ComputeEigenvectors
NumericalIssue
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
ComputationInfo