![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 David Harmon <dharmon@gmail.com> 00005 // 00006 // Eigen is free software; you can redistribute it and/or 00007 // modify it under the terms of the GNU Lesser General Public 00008 // License as published by the Free Software Foundation; either 00009 // version 3 of the License, or (at your option) any later version. 00010 // 00011 // Alternatively, you can redistribute it and/or 00012 // modify it under the terms of the GNU General Public License as 00013 // published by the Free Software Foundation; either version 2 of 00014 // the License, or (at your option) any later version. 00015 // 00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 00019 // GNU General Public License for more details. 00020 // 00021 // You should have received a copy of the GNU Lesser General Public 00022 // License and a copy of the GNU General Public License along with 00023 // Eigen. If not, see <http://www.gnu.org/licenses/>. 00024 00025 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 00026 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 00027 00028 #include <Eigen/Dense> 00029 00030 namespace Eigen { 00031 00032 namespace internal { 00033 template<typename Scalar, typename RealScalar> struct arpack_wrapper; 00034 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> struct OP; 00035 } 00036 00037 00038 00039 template<typename MatrixType, typename MatrixSolver=SimplicialLLT<MatrixType>, bool BisSPD=false> 00040 class ArpackGeneralizedSelfAdjointEigenSolver 00041 { 00042 public: 00043 //typedef typename MatrixSolver::MatrixType MatrixType; 00044 00046 typedef typename MatrixType::Scalar Scalar; 00047 typedef typename MatrixType::Index Index; 00048 00055 typedef typename NumTraits<Scalar>::Real RealScalar; 00056 00062 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; 00063 00070 ArpackGeneralizedSelfAdjointEigenSolver() 00071 : m_eivec(), 00072 m_eivalues(), 00073 m_isInitialized(false), 00074 m_eigenvectorsOk(false), 00075 m_nbrConverged(0), 00076 m_nbrIterations(0) 00077 { } 00078 00101 ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B, 00102 Index nbrEigenvalues, std::string eigs_sigma="LM", 00103 int options=ComputeEigenvectors, RealScalar tol=0.0) 00104 : m_eivec(), 00105 m_eivalues(), 00106 m_isInitialized(false), 00107 m_eigenvectorsOk(false), 00108 m_nbrConverged(0), 00109 m_nbrIterations(0) 00110 { 00111 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 00112 } 00113 00136 ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, 00137 Index nbrEigenvalues, std::string eigs_sigma="LM", 00138 int options=ComputeEigenvectors, RealScalar tol=0.0) 00139 : m_eivec(), 00140 m_eivalues(), 00141 m_isInitialized(false), 00142 m_eigenvectorsOk(false), 00143 m_nbrConverged(0), 00144 m_nbrIterations(0) 00145 { 00146 compute(A, nbrEigenvalues, eigs_sigma, options, tol); 00147 } 00148 00149 00173 ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, const MatrixType& B, 00174 Index nbrEigenvalues, std::string eigs_sigma="LM", 00175 int options=ComputeEigenvectors, RealScalar tol=0.0); 00176 00199 ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, 00200 Index nbrEigenvalues, std::string eigs_sigma="LM", 00201 int options=ComputeEigenvectors, RealScalar tol=0.0); 00202 00203 00223 const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const 00224 { 00225 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 00226 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00227 return m_eivec; 00228 } 00229 00245 const Matrix<Scalar, Dynamic, 1>& eigenvalues() const 00246 { 00247 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 00248 return m_eivalues; 00249 } 00250 00269 Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const 00270 { 00271 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 00272 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00273 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 00274 } 00275 00294 Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const 00295 { 00296 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 00297 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00298 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 00299 } 00300 00305 ComputationInfo info() const 00306 { 00307 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 00308 return m_info; 00309 } 00310 00311 size_t getNbrConvergedEigenValues() const 00312 { return m_nbrConverged; } 00313 00314 size_t getNbrIterations() const 00315 { return m_nbrIterations; } 00316 00317 protected: 00318 Matrix<Scalar, Dynamic, Dynamic> m_eivec; 00319 Matrix<Scalar, Dynamic, 1> m_eivalues; 00320 ComputationInfo m_info; 00321 bool m_isInitialized; 00322 bool m_eigenvectorsOk; 00323 00324 size_t m_nbrConverged; 00325 size_t m_nbrIterations; 00326 }; 00327 00328 00329 00330 00331 00332 template<typename MatrixType, typename MatrixSolver, bool BisSPD> 00333 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 00334 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 00335 ::compute(const MatrixType& A, Index nbrEigenvalues, 00336 std::string eigs_sigma, int options, RealScalar tol) 00337 { 00338 MatrixType B(0,0); 00339 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 00340 00341 return *this; 00342 } 00343 00344 00345 template<typename MatrixType, typename MatrixSolver, bool BisSPD> 00346 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 00347 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 00348 ::compute(const MatrixType& A, const MatrixType& B, Index nbrEigenvalues, 00349 std::string eigs_sigma, int options, RealScalar tol) 00350 { 00351 eigen_assert(A.cols() == A.rows()); 00352 eigen_assert(B.cols() == B.rows()); 00353 eigen_assert(B.rows() == 0 || A.cols() == B.rows()); 00354 eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0 00355 && (options & EigVecMask) != EigVecMask 00356 && "invalid option parameter"); 00357 00358 bool isBempty = (B.rows() == 0) || (B.cols() == 0); 00359 00360 // For clarity, all parameters match their ARPACK name 00361 // 00362 // Always 0 on the first call 00363 // 00364 int ido = 0; 00365 00366 int n = (int)A.cols(); 00367 00368 // User options: "LA", "SA", "SM", "LM", "BE" 00369 // 00370 char whch[3] = "LM"; 00371 00372 // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 } 00373 // 00374 RealScalar sigma = 0.0; 00375 00376 if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) 00377 { 00378 eigs_sigma[0] = toupper(eigs_sigma[0]); 00379 eigs_sigma[1] = toupper(eigs_sigma[1]); 00380 00381 // In the following special case we're going to invert the problem, since solving 00382 // for larger magnitude is much much faster 00383 // i.e., if 'SM' is specified, we're going to really use 'LM', the default 00384 // 00385 if (eigs_sigma.substr(0,2) != "SM") 00386 { 00387 whch[0] = eigs_sigma[0]; 00388 whch[1] = eigs_sigma[1]; 00389 } 00390 } 00391 else 00392 { 00393 eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!"); 00394 00395 // If it's not scalar values, then the user may be explicitly 00396 // specifying the sigma value to cluster the evs around 00397 // 00398 sigma = atof(eigs_sigma.c_str()); 00399 00400 // If atof fails, it returns 0.0, which is a fine default 00401 // 00402 } 00403 00404 // "I" means normal eigenvalue problem, "G" means generalized 00405 // 00406 char bmat[2] = "I"; 00407 if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD)) 00408 bmat[0] = 'G'; 00409 00410 // Now we determine the mode to use 00411 // 00412 int mode = (bmat[0] == 'G') + 1; 00413 if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))) 00414 { 00415 // We're going to use shift-and-invert mode, and basically find 00416 // the largest eigenvalues of the inverse operator 00417 // 00418 mode = 3; 00419 } 00420 00421 // The user-specified number of eigenvalues/vectors to compute 00422 // 00423 int nev = (int)nbrEigenvalues; 00424 00425 // Allocate space for ARPACK to store the residual 00426 // 00427 Scalar *resid = new Scalar[n]; 00428 00429 // Number of Lanczos vectors, must satisfy nev < ncv <= n 00430 // Note that this indicates that nev != n, and we cannot compute 00431 // all eigenvalues of a mtrix 00432 // 00433 int ncv = std::min(std::max(2*nev, 20), n); 00434 00435 // The working n x ncv matrix, also store the final eigenvectors (if computed) 00436 // 00437 Scalar *v = new Scalar[n*ncv]; 00438 int ldv = n; 00439 00440 // Working space 00441 // 00442 Scalar *workd = new Scalar[3*n]; 00443 int lworkl = ncv*ncv+8*ncv; // Must be at least this length 00444 Scalar *workl = new Scalar[lworkl]; 00445 00446 int *iparam= new int[11]; 00447 iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it 00448 iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1))); 00449 iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert 00450 00451 // Used during reverse communicate to notify where arrays start 00452 // 00453 int *ipntr = new int[11]; 00454 00455 // Error codes are returned in here, initial value of 0 indicates a random initial 00456 // residual vector is used, any other values means resid contains the initial residual 00457 // vector, possibly from a previous run 00458 // 00459 int info = 0; 00460 00461 Scalar scale = 1.0; 00462 //if (!isBempty) 00463 //{ 00464 //Scalar scale = B.norm() / std::sqrt(n); 00465 //scale = std::pow(2, std::floor(std::log(scale+1))); 00467 //for (size_t i=0; i<(size_t)B.outerSize(); i++) 00468 // for (typename MatrixType::InnerIterator it(B, i); it; ++it) 00469 // it.valueRef() /= scale; 00470 //} 00471 00472 MatrixSolver OP; 00473 if (mode == 1 || mode == 2) 00474 { 00475 if (!isBempty) 00476 OP.compute(B); 00477 } 00478 else if (mode == 3) 00479 { 00480 if (sigma == 0.0) 00481 { 00482 OP.compute(A); 00483 } 00484 else 00485 { 00486 // Note: We will never enter here because sigma must be 0.0 00487 // 00488 if (isBempty) 00489 { 00490 MatrixType AminusSigmaB(A); 00491 for (Index i=0; i<A.rows(); ++i) 00492 AminusSigmaB.coeffRef(i,i) -= sigma; 00493 00494 OP.compute(AminusSigmaB); 00495 } 00496 else 00497 { 00498 MatrixType AminusSigmaB = A - sigma * B; 00499 OP.compute(AminusSigmaB); 00500 } 00501 } 00502 } 00503 00504 if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success) 00505 std::cout << "Error factoring matrix" << std::endl; 00506 00507 do 00508 { 00509 internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, 00510 &ncv, v, &ldv, iparam, ipntr, workd, workl, 00511 &lworkl, &info); 00512 00513 if (ido == -1 || ido == 1) 00514 { 00515 Scalar *in = workd + ipntr[0] - 1; 00516 Scalar *out = workd + ipntr[1] - 1; 00517 00518 if (ido == 1 && mode != 2) 00519 { 00520 Scalar *out2 = workd + ipntr[2] - 1; 00521 if (isBempty || mode == 1) 00522 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 00523 else 00524 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 00525 00526 in = workd + ipntr[2] - 1; 00527 } 00528 00529 if (mode == 1) 00530 { 00531 if (isBempty) 00532 { 00533 // OP = A 00534 // 00535 Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 00536 } 00537 else 00538 { 00539 // OP = L^{-1}AL^{-T} 00540 // 00541 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out); 00542 } 00543 } 00544 else if (mode == 2) 00545 { 00546 if (ido == 1) 00547 Matrix<Scalar, Dynamic, 1>::Map(in, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 00548 00549 // OP = B^{-1} A 00550 // 00551 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 00552 } 00553 else if (mode == 3) 00554 { 00555 // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I) 00556 // The B * in is already computed and stored at in if ido == 1 00557 // 00558 if (ido == 1 || isBempty) 00559 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 00560 else 00561 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n)); 00562 } 00563 } 00564 else if (ido == 2) 00565 { 00566 Scalar *in = workd + ipntr[0] - 1; 00567 Scalar *out = workd + ipntr[1] - 1; 00568 00569 if (isBempty || mode == 1) 00570 Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 00571 else 00572 Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 00573 } 00574 } while (ido != 99); 00575 00576 if (info == 1) 00577 m_info = NoConvergence; 00578 else if (info == 3) 00579 m_info = NumericalIssue; 00580 else if (info < 0) 00581 m_info = InvalidInput; 00582 else if (info != 0) 00583 eigen_assert(false && "Unknown ARPACK return value!"); 00584 else 00585 { 00586 // Do we compute eigenvectors or not? 00587 // 00588 int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors; 00589 00590 // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK)) 00591 // 00592 char howmny[2] = "A"; 00593 00594 // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK) 00595 // 00596 int *select = new int[ncv]; 00597 00598 // Final eigenvalues 00599 // 00600 m_eivalues.resize(nev, 1); 00601 00602 internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv, 00603 &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv, 00604 v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info); 00605 00606 if (info == -14) 00607 m_info = NoConvergence; 00608 else if (info != 0) 00609 m_info = InvalidInput; 00610 else 00611 { 00612 if (rvec) 00613 { 00614 m_eivec.resize(A.rows(), nev); 00615 for (int i=0; i<nev; i++) 00616 for (int j=0; j<n; j++) 00617 m_eivec(j,i) = v[i*n+j] / scale; 00618 00619 if (mode == 1 && !isBempty && BisSPD) 00620 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data()); 00621 00622 m_eigenvectorsOk = true; 00623 } 00624 00625 m_nbrIterations = iparam[2]; 00626 m_nbrConverged = iparam[4]; 00627 00628 m_info = Success; 00629 } 00630 00631 delete[] select; 00632 } 00633 00634 delete[] v; 00635 delete[] iparam; 00636 delete[] ipntr; 00637 delete[] workd; 00638 delete[] workl; 00639 delete[] resid; 00640 00641 m_isInitialized = true; 00642 00643 return *this; 00644 } 00645 00646 00647 // Single precision 00648 // 00649 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which, 00650 int *nev, float *tol, float *resid, int *ncv, 00651 float *v, int *ldv, int *iparam, int *ipntr, 00652 float *workd, float *workl, int *lworkl, 00653 int *info); 00654 00655 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d, 00656 float *z, int *ldz, float *sigma, 00657 char *bmat, int *n, char *which, int *nev, 00658 float *tol, float *resid, int *ncv, float *v, 00659 int *ldv, int *iparam, int *ipntr, float *workd, 00660 float *workl, int *lworkl, int *ierr); 00661 00662 // Double precision 00663 // 00664 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which, 00665 int *nev, double *tol, double *resid, int *ncv, 00666 double *v, int *ldv, int *iparam, int *ipntr, 00667 double *workd, double *workl, int *lworkl, 00668 int *info); 00669 00670 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d, 00671 double *z, int *ldz, double *sigma, 00672 char *bmat, int *n, char *which, int *nev, 00673 double *tol, double *resid, int *ncv, double *v, 00674 int *ldv, int *iparam, int *ipntr, double *workd, 00675 double *workl, int *lworkl, int *ierr); 00676 00677 00678 namespace internal { 00679 00680 template<typename Scalar, typename RealScalar> struct arpack_wrapper 00681 { 00682 static inline void saupd(int *ido, char *bmat, int *n, char *which, 00683 int *nev, RealScalar *tol, Scalar *resid, int *ncv, 00684 Scalar *v, int *ldv, int *iparam, int *ipntr, 00685 Scalar *workd, Scalar *workl, int *lworkl, int *info) 00686 { 00687 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 00688 } 00689 00690 static inline void seupd(int *rvec, char *All, int *select, Scalar *d, 00691 Scalar *z, int *ldz, RealScalar *sigma, 00692 char *bmat, int *n, char *which, int *nev, 00693 RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, 00694 int *ldv, int *iparam, int *ipntr, Scalar *workd, 00695 Scalar *workl, int *lworkl, int *ierr) 00696 { 00697 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 00698 } 00699 }; 00700 00701 template <> struct arpack_wrapper<float, float> 00702 { 00703 static inline void saupd(int *ido, char *bmat, int *n, char *which, 00704 int *nev, float *tol, float *resid, int *ncv, 00705 float *v, int *ldv, int *iparam, int *ipntr, 00706 float *workd, float *workl, int *lworkl, int *info) 00707 { 00708 ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 00709 } 00710 00711 static inline void seupd(int *rvec, char *All, int *select, float *d, 00712 float *z, int *ldz, float *sigma, 00713 char *bmat, int *n, char *which, int *nev, 00714 float *tol, float *resid, int *ncv, float *v, 00715 int *ldv, int *iparam, int *ipntr, float *workd, 00716 float *workl, int *lworkl, int *ierr) 00717 { 00718 sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 00719 workd, workl, lworkl, ierr); 00720 } 00721 }; 00722 00723 template <> struct arpack_wrapper<double, double> 00724 { 00725 static inline void saupd(int *ido, char *bmat, int *n, char *which, 00726 int *nev, double *tol, double *resid, int *ncv, 00727 double *v, int *ldv, int *iparam, int *ipntr, 00728 double *workd, double *workl, int *lworkl, int *info) 00729 { 00730 dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 00731 } 00732 00733 static inline void seupd(int *rvec, char *All, int *select, double *d, 00734 double *z, int *ldz, double *sigma, 00735 char *bmat, int *n, char *which, int *nev, 00736 double *tol, double *resid, int *ncv, double *v, 00737 int *ldv, int *iparam, int *ipntr, double *workd, 00738 double *workl, int *lworkl, int *ierr) 00739 { 00740 dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 00741 workd, workl, lworkl, ierr); 00742 } 00743 }; 00744 00745 00746 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> 00747 struct OP 00748 { 00749 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out); 00750 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs); 00751 }; 00752 00753 template<typename MatrixSolver, typename MatrixType, typename Scalar> 00754 struct OP<MatrixSolver, MatrixType, Scalar, true> 00755 { 00756 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 00757 { 00758 // OP = L^{-1} A L^{-T} (B = LL^T) 00759 // 00760 // First solve L^T out = in 00761 // 00762 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 00763 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 00764 00765 // Then compute out = A out 00766 // 00767 Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n); 00768 00769 // Then solve L out = out 00770 // 00771 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 00772 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n)); 00773 } 00774 00775 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 00776 { 00777 // Solve L^T out = in 00778 // 00779 Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k)); 00780 Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k); 00781 } 00782 00783 }; 00784 00785 template<typename MatrixSolver, typename MatrixType, typename Scalar> 00786 struct OP<MatrixSolver, MatrixType, Scalar, false> 00787 { 00788 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 00789 { 00790 eigen_assert(false && "Should never be in here..."); 00791 } 00792 00793 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 00794 { 00795 eigen_assert(false && "Should never be in here..."); 00796 } 00797 00798 }; 00799 00800 } // end namespace internal 00801 00802 } // end namespace Eigen 00803 00804 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H 00805