Tapkee
arpack_wrapper.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * This code is based on Eigen3 ARPACK wrapper
00004  * Copyright (c) 2012 David Harmon
00005  *
00006  * The code is relicensed under BSD 3-clause with the permission of its
00007  * original author. Some changes were made afterwards by Sergey Lisitsyn.
00008  */
00009 
00010 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
00011 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H
00012 
00013 /* Tapkee includes */
00014 #include <tapkee/defines.hpp>
00015 /* End of Tapkee includes */
00016 
00017 using Eigen::Matrix;
00018 using Eigen::Dynamic;
00019 using Eigen::ComputationInfo;
00020 using Eigen::ComputeEigenvectors;
00021 using Eigen::NoConvergence;
00022 using Eigen::NumericalIssue;
00023 #ifndef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
00024 using Eigen::InvalidInput;
00025 #endif
00026 using Eigen::Success;
00027 using Eigen::NumTraits;
00028 
00030 namespace arpack_internal 
00031 {
00032     template<typename Scalar, typename RealScalar> struct arpack_wrapper;
00033 }
00034 
00035 template<class LMatrixType, class RMatrixType, class MatrixOperation, bool BisSPD=false>
00036 class ArpackGeneralizedSelfAdjointEigenSolver
00037 {
00038 public:
00039   //typedef typename MatrixSolver::MatrixType MatrixType;
00040 
00042     typedef typename LMatrixType::Scalar Scalar;
00043     typedef typename LMatrixType::Index Index;
00044 
00051     typedef typename NumTraits<Scalar>::Real RealScalar;
00052 
00058     typedef typename Eigen::internal::plain_col_type<LMatrixType, RealScalar>::type RealVectorType;
00059 
00066     ArpackGeneralizedSelfAdjointEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false),
00067         m_eigenvectorsOk(false), m_nbrConverged(0), m_nbrIterations(0)
00068     { 
00069     }
00070 
00093     ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, const RMatrixType& B,
00094                                             Index nbrEigenvalues, std::string eigs_sigma="LM",
00095                                             int parameters=ComputeEigenvectors, RealScalar tol=0.0) : 
00096         m_eivec(), m_eivalues(), m_info(), m_isInitialized(false), m_eigenvectorsOk(false),
00097         m_nbrConverged(0), m_nbrIterations(0)
00098     {
00099         compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
00100     }
00101 
00123     ArpackGeneralizedSelfAdjointEigenSolver(const LMatrixType& A, Index nbrEigenvalues, std::string eigs_sigma="LM",
00124                                             int parameters=ComputeEigenvectors, RealScalar tol=0.0) : 
00125         m_eivec(), m_eivalues(), m_info(), m_isInitialized(false), m_eigenvectorsOk(false),
00126         m_nbrConverged(0), m_nbrIterations(0)
00127     {
00128         compute(A, nbrEigenvalues, eigs_sigma, parameters, tol);
00129     }
00130 
00154     ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, const RMatrixType& B,
00155                                                      Index nbrEigenvalues, std::string eigs_sigma="LM",
00156                                                      int parameters=ComputeEigenvectors, RealScalar tol=0.0);
00157 
00180     ArpackGeneralizedSelfAdjointEigenSolver& compute(const LMatrixType& A, Index nbrEigenvalues, 
00181                                                      std::string eigs_sigma="LM",
00182                                                      int parameters=ComputeEigenvectors, RealScalar tol=0.0);
00183 
00200     const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const
00201     {
00202         eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
00203         eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00204         return m_eivec;
00205     }
00206 
00219     const Matrix<Scalar, Dynamic, 1>& eigenvalues() const
00220     {
00221         eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
00222         return m_eivalues;
00223     }
00224 
00239     Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const
00240     {
00241         eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00242         eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00243         return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00244     }
00245 
00260     Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const
00261     {
00262         eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00263         eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00264         return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00265     }
00266 
00271     ComputationInfo info() const
00272     {
00273         eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized.");
00274         return m_info;
00275     }
00276 
00277     size_t getNbrConvergedEigenValues() const
00278     {
00279         return m_nbrConverged;
00280     }
00281 
00282     size_t getNbrIterations() const
00283     {
00284         return m_nbrIterations; 
00285     }
00286 
00287 protected:
00288 
00289     Matrix<Scalar, Dynamic, Dynamic> m_eivec;
00290     Matrix<Scalar, Dynamic, 1> m_eivalues;
00291     ComputationInfo m_info;
00292     bool m_isInitialized;
00293     bool m_eigenvectorsOk;
00294 
00295     size_t m_nbrConverged;
00296     size_t m_nbrIterations;
00297 };
00298 
00299 template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
00300 ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>&
00301     ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>
00302 ::compute(const LMatrixType& A, Index nbrEigenvalues,
00303           std::string eigs_sigma, int parameters, RealScalar tol)
00304 {
00305     RMatrixType B(0,0);
00306     return compute(A, B, nbrEigenvalues, eigs_sigma, parameters, tol);
00307 }
00308 
00309 template<typename LMatrixType, typename RMatrixType, class MatrixOperation, bool BisSPD>
00310 ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>&
00311     ArpackGeneralizedSelfAdjointEigenSolver<LMatrixType, RMatrixType, MatrixOperation, BisSPD>
00312 ::compute(const LMatrixType& A, const RMatrixType& B, Index nbrEigenvalues,
00313           std::string eigs_sigma, int parameters, RealScalar tol)
00314 {
00315     eigen_assert(B.cols() == B.rows());
00316     eigen_assert(B.rows() == 0 || A.rows() == B.rows() || A.cols() == B.cols());
00317     eigen_assert((parameters &~ (Eigen::EigVecMask | Eigen::GenEigMask)) == 0
00318                  && (parameters & Eigen::EigVecMask) != Eigen::EigVecMask
00319                  && "invalid option parameter");
00320 
00321     bool isBempty = (B.rows() == 0) || (B.cols() == 0);
00322 
00323     // For clarity, all parameters match their ARPACK name
00324     // Always 0 on the first call
00325     int ido = 0;
00326     int n = (int)A.rows();
00327     // User parameters: "LA", "SA", "SM", "LM", "BE"
00328     char whch[3] = "LM";
00329 
00330     // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 }
00331     RealScalar sigma = 0.0;
00332 
00333     if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))
00334     {
00335         eigs_sigma[0] = toupper(eigs_sigma[0]);
00336         eigs_sigma[1] = toupper(eigs_sigma[1]);
00337         // In the following special case we're going to invert the problem, since solving
00338         // for larger magnitude is much much faster
00339         // i.e., if 'SM' is specified, we're going to really use 'LM', the default
00340         if (eigs_sigma.substr(0,2) != "SM")
00341         {
00342             whch[0] = eigs_sigma[0];
00343             whch[1] = eigs_sigma[1];
00344         }
00345     }
00346     else
00347     {
00348         eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!");
00349         // If it's not scalar values, then the user may be explicitly
00350         // specifying the sigma value to cluster the evs around
00351         sigma = atof(eigs_sigma.c_str());
00352         // If atof fails, it returns 0.0, which is a fine default
00353     }
00354 
00355     // "I" means normal eigenvalue problem, "G" means generalized
00356     //
00357     char bmat[2] = "I";
00358     if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD))
00359         bmat[0] = 'G';
00360 
00361     // Now we determine the mode to use
00362     //
00363     int mode = (bmat[0] == 'G') + 1;
00364     if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])))
00365     {
00366         // We're going to use shift-and-invert mode, and basically find
00367         // the largest eigenvalues of the inverse operator
00368         mode = 3;
00369         //std::cout << "Mode = 3\n";
00370     }
00371 
00372     // The user-specified number of eigenvalues/vectors to compute
00373     int nev = (int)nbrEigenvalues;
00374     // Allocate space for ARPACK to store the residual
00375     Scalar *resid = new Scalar[n];
00376     // Number of Lanczos vectors, must satisfy nev < ncv <= n
00377     // Note that this indicates that nev != n, and we cannot compute
00378     // all eigenvalues of a mtrix
00379     int ncv = std::min(std::max(2*nev, 20), n);
00380 
00381     // The working n x ncv matrix, also store the final eigenvectors (if computed)
00382     Scalar *v = new Scalar[n*ncv];
00383     int ldv = n;
00384 
00385     // Working space
00386     Scalar *workd = new Scalar[3*n];
00387     int lworkl = ncv*ncv+8*ncv; // Must be at least this length
00388     Scalar *workl = new Scalar[lworkl];
00389 
00390     int *iparam= new int[11];
00391     iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it
00392     iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1)));
00393     iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert
00394 
00395     // Used during reverse communicate to notify where arrays start
00396     int *ipntr = new int[11]; 
00397 
00398     // Error codes are returned in here, initial value of 0 indicates a random initial
00399     // residual vector is used, any other values means resid contains the initial residual
00400     // vector, possibly from a previous run
00401     int cinfo = 0;
00402 
00403     Scalar scale = 1.0;
00404     /*
00405     if (!isBempty)
00406     {
00407     Scalar scale = B.norm() / std::sqrt(n);
00408     scale = std::pow(2, std::floor(std::log(scale+1)));
00410     for (size_t i=0; i<(size_t)B.outerSize(); i++)
00411         for (typename RMatrixType::InnerIterator it(B, i); it; ++it)
00412             it.valueRef() /= scale;
00413     }
00414     */
00415 
00416     MatrixOperation op(A);
00417     //      (mode==1 || mode==2) ? B :
00418     //                   (mode==3) ? A : LMatrixType());
00419 
00420     do
00421     {
00422         //std::cout << "Entering main loop\n";
00423         arpack_internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, 
00424                                                                    &ncv, v, &ldv, iparam, ipntr, workd, workl,
00425                                                                    &lworkl, &cinfo);
00426         if (ido == -1 || ido == 1)
00427         {
00428             Scalar *in  = workd + ipntr[0] - 1;
00429             Scalar *out = workd + ipntr[1] - 1;
00430 
00431             if (ido == 1 && mode != 2)
00432             {
00433                 Scalar *out2 = workd + ipntr[2] - 1;
00434                 if (isBempty || mode == 1)
00435                     Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
00436                 else
00437                     Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
00438                 in = workd + ipntr[2] - 1;
00439             }
00440 
00441             if (mode == 1)
00442             {
00443                 if (isBempty)
00444                 {
00445                     // OP = A
00446                     Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
00447                 }
00448                 else
00449                 {
00450                     // TODO OP = L^{-1}AL^{-T}
00451                     //internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out);
00452                 }
00453             }
00454             else if (mode == 2)
00455             {
00456                 if (ido == 1)
00457                     Matrix<Scalar, Dynamic, 1>::Map(in, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
00458                 // OP = B^{-1} A
00459                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
00460             }
00461             else if (mode == 3)
00462             {
00463                 // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I)
00464                 // The B * in is already computed and stored at in if ido == 1
00465                 if (ido == 1 || isBempty)
00466                     Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(Matrix<Scalar, Dynamic, 1>::Map(in, n));
00467                 else
00468                     Matrix<Scalar, Dynamic, 1>::Map(out, n) = op(B * Matrix<Scalar, Dynamic, 1>::Map(in, n));
00469             }
00470         }
00471         else if (ido == 2)
00472         {
00473             Scalar *in  = workd + ipntr[0] - 1;
00474             Scalar *out = workd + ipntr[1] - 1;
00475 
00476             if (isBempty || mode == 1)
00477                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n);
00478             else
00479                 Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n);
00480         }
00481         //Scalar *out = workd + ipntr[1] - 1;
00482         //std::cout << Matrix<Scalar, Dynamic, 1>::Map(out, n).transpose() << std::endl;
00483     } while (ido != 99);
00484         
00485     //std::cout << "Finsihed\n";
00486 
00487     if (cinfo == 1)
00488     {
00489         //std::cout << "FAILED WITH NO CONV\n";
00490         m_info = NoConvergence;
00491     }
00492     else if (cinfo == 3)
00493     {
00494         //std::cout << "FAILED WITH NUMERICAL ISSUE\n";
00495         m_info = NumericalIssue;
00496     }
00497     else if (cinfo < 0)
00498     {
00499         //std::cout << "FAILED WITH INVALID INPUT " << info << "\n";
00500 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
00501         m_info = NoConvergence;
00502 #else
00503         m_info = InvalidInput;
00504 #endif
00505     }
00506     else if (cinfo != 0)
00507         eigen_assert(false && "Unknown ARPACK return value!");
00508     else
00509     {
00510         // Do we compute eigenvectors or not?
00511         int rvec = (parameters & ComputeEigenvectors) == ComputeEigenvectors;
00512 
00513         // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK))
00514         char howmny[2] = "A"; 
00515 
00516         // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK)
00517         int *select = new int[ncv];
00518 
00519         // Final eigenvalues
00520         m_eivalues.resize(nev, 1);
00521 
00522         arpack_internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv,
00523                                                                    &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv,
00524                                                                    v, &ldv, iparam, ipntr, workd, workl, &lworkl, &cinfo);
00525 
00526         if (cinfo == -14)
00527             m_info = NoConvergence;
00528         else if (cinfo != 0)
00529 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
00530             m_info = NoConvergence;
00531 #else
00532             m_info = InvalidInput;
00533 #endif
00534         else
00535         {
00536             if (rvec)
00537             {
00538                 m_eivec.resize(A.rows(), nev);
00539                 for (int i=0; i<nev; i++)
00540                     for (int j=0; j<n; j++)
00541                         m_eivec(j,i) = v[i*n+j] / scale;
00542   
00543                 // TODO if (mode == 1 && !isBempty && BisSPD)
00544                 //  internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data());
00545 
00546                 m_eigenvectorsOk = true;
00547             }
00548 
00549             m_nbrIterations = iparam[2];
00550             m_nbrConverged  = iparam[4];
00551             m_info = Success;
00552         }
00553         delete[] select;
00554     }
00555 
00556     delete[] v;
00557     delete[] iparam;
00558     delete[] ipntr;
00559     delete[] workd;
00560     delete[] workl;
00561     delete[] resid;
00562     m_isInitialized = true;
00563     return *this;
00564 }
00565 
00566 // Single precision
00567 //
00568 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which,
00569                         int *nev, float *tol, float *resid, int *ncv,
00570                         float *v, int *ldv, int *iparam, int *ipntr,
00571                         float *workd, float *workl, int *lworkl,
00572                         int *info);
00573 
00574 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d,
00575                         float *z, int *ldz, float *sigma, 
00576                         char *bmat, int *n, char *which, int *nev,
00577                         float *tol, float *resid, int *ncv, float *v,
00578                         int *ldv, int *iparam, int *ipntr, float *workd,
00579                         float *workl, int *lworkl, int *ierr);
00580 
00581 // Double precision
00582 //
00583 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which,
00584                         int *nev, double *tol, double *resid, int *ncv,
00585                         double *v, int *ldv, int *iparam, int *ipntr,
00586                         double *workd, double *workl, int *lworkl,
00587                         int *info);
00588 
00589 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d,
00590                         double *z, int *ldz, double *sigma, 
00591                         char *bmat, int *n, char *which, int *nev,
00592                         double *tol, double *resid, int *ncv, double *v,
00593                         int *ldv, int *iparam, int *ipntr, double *workd,
00594                         double *workl, int *lworkl, int *ierr);
00595 
00596 namespace arpack_internal {
00597 
00598 template<typename Scalar, typename RealScalar> struct arpack_wrapper
00599 {
00600   static inline void saupd(int *ido, char *bmat, int *n, char *which,
00601       int *nev, RealScalar *tol, Scalar *resid, int *ncv,
00602       Scalar *v, int *ldv, int *iparam, int *ipntr,
00603       Scalar *workd, Scalar *workl, int *lworkl, int *info);
00604 
00605   static inline void seupd(int *rvec, char *All, int *select, Scalar *d,
00606       Scalar *z, int *ldz, RealScalar *sigma,
00607       char *bmat, int *n, char *which, int *nev,
00608       RealScalar *tol, Scalar *resid, int *ncv, Scalar *v,
00609       int *ldv, int *iparam, int *ipntr, Scalar *workd,
00610       Scalar *workl, int *lworkl, int *ierr);
00611 };
00612 
00613 template <> struct arpack_wrapper<float, float>
00614 {
00615   static inline void saupd(int *ido, char *bmat, int *n, char *which,
00616       int *nev, float *tol, float *resid, int *ncv,
00617       float *v, int *ldv, int *iparam, int *ipntr,
00618       float *workd, float *workl, int *lworkl, int *info)
00619   {
00620     ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
00621   }
00622 
00623   static inline void seupd(int *rvec, char *All, int *select, float *d,
00624       float* z, int* ldz, float *sigma,
00625       char *bmat, int *n, char *which, int *nev,
00626       float *tol, float *resid, int *ncv, float *v,
00627       int *ldv, int *iparam, int *ipntr, float *workd,
00628       float *workl, int *lworkl, int *ierr)
00629   {
00630     sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
00631         workd, workl, lworkl, ierr);
00632   }
00633 };
00634 
00635 template <> struct arpack_wrapper<double, double>
00636 {
00637   static inline void saupd(int *ido, char *bmat, int *n, char *which,
00638       int *nev, double *tol, double *resid, int *ncv,
00639       double *v, int *ldv, int *iparam, int *ipntr,
00640       double *workd, double *workl, int *lworkl, int *info)
00641   {
00642     dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
00643   }
00644 
00645   static inline void seupd(int *rvec, char *All, int *select, double *d,
00646       double* z, int* ldz, double *sigma,
00647       char *bmat, int *n, char *which, int *nev,
00648       double *tol, double *resid, int *ncv, double *v,
00649       int *ldv, int *iparam, int *ipntr, double *workd,
00650       double *workl, int *lworkl, int *ierr)
00651   {
00652     dseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr,
00653         workd, workl, lworkl, ierr);
00654   }
00655 };
00656 
00657 }
00658 
00659 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H
00660 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines