ArpackSelfAdjointEigenSolver.h (29075B)
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 30 /** \brief Scalar type for matrices of type \p MatrixType. */ 31 typedef typename MatrixType::Scalar Scalar; 32 typedef typename MatrixType::Index Index; 33 34 /** \brief Real scalar type for \p MatrixType. 35 * 36 * This is just \c Scalar if #Scalar is real (e.g., \c float or 37 * \c Scalar), and the type of the real part of \c Scalar if #Scalar is 38 * complex. 39 */ 40 typedef typename NumTraits<Scalar>::Real RealScalar; 41 42 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 43 * 44 * This is a column vector with entries of type #RealScalar. 45 * The length of the vector is the size of \p nbrEigenvalues. 46 */ 47 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; 48 49 /** \brief Default constructor. 50 * 51 * The default constructor is for cases in which the user intends to 52 * perform decompositions via compute(). 53 * 54 */ 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 64 /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix. 65 * 66 * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 67 * computed. By default, the upper triangular part is used, but can be changed 68 * through the template parameter. 69 * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem. 70 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 71 * Must be less than the size of the input matrix, or an error is returned. 72 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 73 * respective meanings to find the largest magnitude , smallest magnitude, 74 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 75 * value can contain floating point value in string form, in which case the 76 * eigenvalues closest to this value will be found. 77 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 78 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 79 * means machine precision. 80 * 81 * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar) 82 * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if 83 * \p options equals #ComputeEigenvectors. 84 * 85 */ 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 99 /** \brief Constructor; computes eigenvalues of given matrix. 100 * 101 * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 102 * computed. By default, the upper triangular part is used, but can be changed 103 * through the template parameter. 104 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 105 * Must be less than the size of the input matrix, or an error is returned. 106 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 107 * respective meanings to find the largest magnitude , smallest magnitude, 108 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 109 * value can contain floating point value in string form, in which case the 110 * eigenvalues closest to this value will be found. 111 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 112 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 113 * means machine precision. 114 * 115 * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar) 116 * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if 117 * \p options equals #ComputeEigenvectors. 118 * 119 */ 120 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 135 /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library. 136 * 137 * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 138 * \param[in] B Selfadjoint matrix for generalized eigenvalues. 139 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 140 * Must be less than the size of the input matrix, or an error is returned. 141 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 142 * respective meanings to find the largest magnitude , smallest magnitude, 143 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 144 * value can contain floating point value in string form, in which case the 145 * eigenvalues closest to this value will be found. 146 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 147 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 148 * means machine precision. 149 * 150 * \returns Reference to \c *this 151 * 152 * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK. The eigenvalues() 153 * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 154 * then the eigenvectors are also computed and can be retrieved by 155 * calling eigenvectors(). 156 * 157 */ 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 162 /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library. 163 * 164 * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 165 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 166 * Must be less than the size of the input matrix, or an error is returned. 167 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 168 * respective meanings to find the largest magnitude , smallest magnitude, 169 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 170 * value can contain floating point value in string form, in which case the 171 * eigenvalues closest to this value will be found. 172 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 173 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 174 * means machine precision. 175 * 176 * \returns Reference to \c *this 177 * 178 * This function computes the eigenvalues of \p A using ARPACK. The eigenvalues() 179 * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 180 * then the eigenvectors are also computed and can be retrieved by 181 * calling eigenvectors(). 182 * 183 */ 184 ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, 185 Index nbrEigenvalues, std::string eigs_sigma="LM", 186 int options=ComputeEigenvectors, RealScalar tol=0.0); 187 188 189 /** \brief Returns the eigenvectors of given matrix. 190 * 191 * \returns A const reference to the matrix whose columns are the eigenvectors. 192 * 193 * \pre The eigenvectors have been computed before. 194 * 195 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 196 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 197 * eigenvectors are normalized to have (Euclidean) norm equal to one. If 198 * this object was used to solve the eigenproblem for the selfadjoint 199 * matrix \f$ A \f$, then the matrix returned by this function is the 200 * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$. 201 * For the generalized eigenproblem, the matrix returned is the solution \f$ A V = D B V \f$ 202 * 203 * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp 204 * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out 205 * 206 * \sa eigenvalues() 207 */ 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 215 /** \brief Returns the eigenvalues of given matrix. 216 * 217 * \returns A const reference to the column vector containing the eigenvalues. 218 * 219 * \pre The eigenvalues have been computed before. 220 * 221 * The eigenvalues are repeated according to their algebraic multiplicity, 222 * so there are as many eigenvalues as rows in the matrix. The eigenvalues 223 * are sorted in increasing order. 224 * 225 * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp 226 * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out 227 * 228 * \sa eigenvectors(), MatrixBase::eigenvalues() 229 */ 230 const Matrix<Scalar, Dynamic, 1>& eigenvalues() const 231 { 232 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 233 return m_eivalues; 234 } 235 236 /** \brief Computes the positive-definite square root of the matrix. 237 * 238 * \returns the positive-definite square root of the matrix 239 * 240 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 241 * have been computed before. 242 * 243 * The square root of a positive-definite matrix \f$ A \f$ is the 244 * positive-definite matrix whose square equals \f$ A \f$. This function 245 * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the 246 * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$. 247 * 248 * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp 249 * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out 250 * 251 * \sa operatorInverseSqrt(), 252 * \ref MatrixFunctions_Module "MatrixFunctions Module" 253 */ 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 261 /** \brief Computes the inverse square root of the matrix. 262 * 263 * \returns the inverse positive-definite square root of the matrix 264 * 265 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 266 * have been computed before. 267 * 268 * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to 269 * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is 270 * cheaper than first computing the square root with operatorSqrt() and 271 * then its inverse with MatrixBase::inverse(). 272 * 273 * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp 274 * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out 275 * 276 * \sa operatorSqrt(), MatrixBase::inverse(), 277 * \ref MatrixFunctions_Module "MatrixFunctions Module" 278 */ 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 286 /** \brief Reports whether previous computation was successful. 287 * 288 * \returns \c Success if computation was successful, \c NoConvergence otherwise. 289 */ 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))); 451 ////M /= scale; 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