1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 David Harmon <dharmon (at) gmail.com> 5 // 6 // Eigen is free software; you can redistribute it and/or 7 // modify it under the terms of the GNU Lesser General Public 8 // License as published by the Free Software Foundation; either 9 // version 3 of the License, or (at your option) any later version. 10 // 11 // Alternatively, you can redistribute it and/or 12 // modify it under the terms of the GNU General Public License as 13 // published by the Free Software Foundation; either version 2 of 14 // the License, or (at your option) any later version. 15 // 16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU Lesser General Public 22 // License and a copy of the GNU General Public License along with 23 // Eigen. If not, see <http://www.gnu.org/licenses/>. 24 25 #ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 26 #define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 27 28 #include <Eigen/Dense> 29 30 namespace Eigen { 31 32 namespace internal { 33 template<typename Scalar, typename RealScalar> struct arpack_wrapper; 34 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> struct OP; 35 } 36 37 38 39 template<typename MatrixType, typename MatrixSolver=SimplicialLLT<MatrixType>, bool BisSPD=false> 40 class ArpackGeneralizedSelfAdjointEigenSolver 41 { 42 public: 43 //typedef typename MatrixSolver::MatrixType MatrixType; 44 45 /** \brief Scalar type for matrices of type \p MatrixType. */ 46 typedef typename MatrixType::Scalar Scalar; 47 typedef typename MatrixType::Index Index; 48 49 /** \brief Real scalar type for \p MatrixType. 50 * 51 * This is just \c Scalar if #Scalar is real (e.g., \c float or 52 * \c Scalar), and the type of the real part of \c Scalar if #Scalar is 53 * complex. 54 */ 55 typedef typename NumTraits<Scalar>::Real RealScalar; 56 57 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 58 * 59 * This is a column vector with entries of type #RealScalar. 60 * The length of the vector is the size of \p nbrEigenvalues. 61 */ 62 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; 63 64 /** \brief Default constructor. 65 * 66 * The default constructor is for cases in which the user intends to 67 * perform decompositions via compute(). 68 * 69 */ 70 ArpackGeneralizedSelfAdjointEigenSolver() 71 : m_eivec(), 72 m_eivalues(), 73 m_isInitialized(false), 74 m_eigenvectorsOk(false), 75 m_nbrConverged(0), 76 m_nbrIterations(0) 77 { } 78 79 /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix. 80 * 81 * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 82 * computed. By default, the upper triangular part is used, but can be changed 83 * through the template parameter. 84 * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem. 85 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 86 * Must be less than the size of the input matrix, or an error is returned. 87 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 88 * respective meanings to find the largest magnitude , smallest magnitude, 89 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 90 * value can contain floating point value in string form, in which case the 91 * eigenvalues closest to this value will be found. 92 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 93 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 94 * means machine precision. 95 * 96 * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar) 97 * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if 98 * \p options equals #ComputeEigenvectors. 99 * 100 */ 101 ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B, 102 Index nbrEigenvalues, std::string eigs_sigma="LM", 103 int options=ComputeEigenvectors, RealScalar tol=0.0) 104 : m_eivec(), 105 m_eivalues(), 106 m_isInitialized(false), 107 m_eigenvectorsOk(false), 108 m_nbrConverged(0), 109 m_nbrIterations(0) 110 { 111 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 112 } 113 114 /** \brief Constructor; computes eigenvalues of given matrix. 115 * 116 * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 117 * computed. By default, the upper triangular part is used, but can be changed 118 * through the template parameter. 119 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 120 * Must be less than the size of the input matrix, or an error is returned. 121 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 122 * respective meanings to find the largest magnitude , smallest magnitude, 123 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 124 * value can contain floating point value in string form, in which case the 125 * eigenvalues closest to this value will be found. 126 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 127 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 128 * means machine precision. 129 * 130 * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar) 131 * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if 132 * \p options equals #ComputeEigenvectors. 133 * 134 */ 135 136 ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, 137 Index nbrEigenvalues, std::string eigs_sigma="LM", 138 int options=ComputeEigenvectors, RealScalar tol=0.0) 139 : m_eivec(), 140 m_eivalues(), 141 m_isInitialized(false), 142 m_eigenvectorsOk(false), 143 m_nbrConverged(0), 144 m_nbrIterations(0) 145 { 146 compute(A, nbrEigenvalues, eigs_sigma, options, tol); 147 } 148 149 150 /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library. 151 * 152 * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 153 * \param[in] B Selfadjoint matrix for generalized eigenvalues. 154 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 155 * Must be less than the size of the input matrix, or an error is returned. 156 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 157 * respective meanings to find the largest magnitude , smallest magnitude, 158 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 159 * value can contain floating point value in string form, in which case the 160 * eigenvalues closest to this value will be found. 161 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 162 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 163 * means machine precision. 164 * 165 * \returns Reference to \c *this 166 * 167 * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK. The eigenvalues() 168 * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 169 * then the eigenvectors are also computed and can be retrieved by 170 * calling eigenvectors(). 171 * 172 */ 173 ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, const MatrixType& B, 174 Index nbrEigenvalues, std::string eigs_sigma="LM", 175 int options=ComputeEigenvectors, RealScalar tol=0.0); 176 177 /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library. 178 * 179 * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 180 * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 181 * Must be less than the size of the input matrix, or an error is returned. 182 * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 183 * respective meanings to find the largest magnitude , smallest magnitude, 184 * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 185 * value can contain floating point value in string form, in which case the 186 * eigenvalues closest to this value will be found. 187 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 188 * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 189 * means machine precision. 190 * 191 * \returns Reference to \c *this 192 * 193 * This function computes the eigenvalues of \p A using ARPACK. The eigenvalues() 194 * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 195 * then the eigenvectors are also computed and can be retrieved by 196 * calling eigenvectors(). 197 * 198 */ 199 ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, 200 Index nbrEigenvalues, std::string eigs_sigma="LM", 201 int options=ComputeEigenvectors, RealScalar tol=0.0); 202 203 204 /** \brief Returns the eigenvectors of given matrix. 205 * 206 * \returns A const reference to the matrix whose columns are the eigenvectors. 207 * 208 * \pre The eigenvectors have been computed before. 209 * 210 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 211 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 212 * eigenvectors are normalized to have (Euclidean) norm equal to one. If 213 * this object was used to solve the eigenproblem for the selfadjoint 214 * matrix \f$ A \f$, then the matrix returned by this function is the 215 * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$. 216 * For the generalized eigenproblem, the matrix returned is the solution \f$ A V = D B V \f$ 217 * 218 * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp 219 * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out 220 * 221 * \sa eigenvalues() 222 */ 223 const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const 224 { 225 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 226 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 227 return m_eivec; 228 } 229 230 /** \brief Returns the eigenvalues of given matrix. 231 * 232 * \returns A const reference to the column vector containing the eigenvalues. 233 * 234 * \pre The eigenvalues have been computed before. 235 * 236 * The eigenvalues are repeated according to their algebraic multiplicity, 237 * so there are as many eigenvalues as rows in the matrix. The eigenvalues 238 * are sorted in increasing order. 239 * 240 * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp 241 * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out 242 * 243 * \sa eigenvectors(), MatrixBase::eigenvalues() 244 */ 245 const Matrix<Scalar, Dynamic, 1>& eigenvalues() const 246 { 247 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 248 return m_eivalues; 249 } 250 251 /** \brief Computes the positive-definite square root of the matrix. 252 * 253 * \returns the positive-definite square root of the matrix 254 * 255 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 256 * have been computed before. 257 * 258 * The square root of a positive-definite matrix \f$ A \f$ is the 259 * positive-definite matrix whose square equals \f$ A \f$. This function 260 * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the 261 * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$. 262 * 263 * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp 264 * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out 265 * 266 * \sa operatorInverseSqrt(), 267 * \ref MatrixFunctions_Module "MatrixFunctions Module" 268 */ 269 Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const 270 { 271 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 272 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 273 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 274 } 275 276 /** \brief Computes the inverse square root of the matrix. 277 * 278 * \returns the inverse positive-definite square root of the matrix 279 * 280 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 281 * have been computed before. 282 * 283 * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to 284 * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is 285 * cheaper than first computing the square root with operatorSqrt() and 286 * then its inverse with MatrixBase::inverse(). 287 * 288 * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp 289 * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out 290 * 291 * \sa operatorSqrt(), MatrixBase::inverse(), 292 * \ref MatrixFunctions_Module "MatrixFunctions Module" 293 */ 294 Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const 295 { 296 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 297 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 298 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 299 } 300 301 /** \brief Reports whether previous computation was successful. 302 * 303 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 304 */ 305 ComputationInfo info() const 306 { 307 eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 308 return m_info; 309 } 310 311 size_t getNbrConvergedEigenValues() const 312 { return m_nbrConverged; } 313 314 size_t getNbrIterations() const 315 { return m_nbrIterations; } 316 317 protected: 318 Matrix<Scalar, Dynamic, Dynamic> m_eivec; 319 Matrix<Scalar, Dynamic, 1> m_eivalues; 320 ComputationInfo m_info; 321 bool m_isInitialized; 322 bool m_eigenvectorsOk; 323 324 size_t m_nbrConverged; 325 size_t m_nbrIterations; 326 }; 327 328 329 330 331 332 template<typename MatrixType, typename MatrixSolver, bool BisSPD> 333 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 334 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 335 ::compute(const MatrixType& A, Index nbrEigenvalues, 336 std::string eigs_sigma, int options, RealScalar tol) 337 { 338 MatrixType B(0,0); 339 compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 340 341 return *this; 342 } 343 344 345 template<typename MatrixType, typename MatrixSolver, bool BisSPD> 346 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 347 ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 348 ::compute(const MatrixType& A, const MatrixType& B, Index nbrEigenvalues, 349 std::string eigs_sigma, int options, RealScalar tol) 350 { 351 eigen_assert(A.cols() == A.rows()); 352 eigen_assert(B.cols() == B.rows()); 353 eigen_assert(B.rows() == 0 || A.cols() == B.rows()); 354 eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0 355 && (options & EigVecMask) != EigVecMask 356 && "invalid option parameter"); 357 358 bool isBempty = (B.rows() == 0) || (B.cols() == 0); 359 360 // For clarity, all parameters match their ARPACK name 361 // 362 // Always 0 on the first call 363 // 364 int ido = 0; 365 366 int n = (int)A.cols(); 367 368 // User options: "LA", "SA", "SM", "LM", "BE" 369 // 370 char whch[3] = "LM"; 371 372 // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 } 373 // 374 RealScalar sigma = 0.0; 375 376 if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) 377 { 378 eigs_sigma[0] = toupper(eigs_sigma[0]); 379 eigs_sigma[1] = toupper(eigs_sigma[1]); 380 381 // In the following special case we're going to invert the problem, since solving 382 // for larger magnitude is much much faster 383 // i.e., if 'SM' is specified, we're going to really use 'LM', the default 384 // 385 if (eigs_sigma.substr(0,2) != "SM") 386 { 387 whch[0] = eigs_sigma[0]; 388 whch[1] = eigs_sigma[1]; 389 } 390 } 391 else 392 { 393 eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!"); 394 395 // If it's not scalar values, then the user may be explicitly 396 // specifying the sigma value to cluster the evs around 397 // 398 sigma = atof(eigs_sigma.c_str()); 399 400 // If atof fails, it returns 0.0, which is a fine default 401 // 402 } 403 404 // "I" means normal eigenvalue problem, "G" means generalized 405 // 406 char bmat[2] = "I"; 407 if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD)) 408 bmat[0] = 'G'; 409 410 // Now we determine the mode to use 411 // 412 int mode = (bmat[0] == 'G') + 1; 413 if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))) 414 { 415 // We're going to use shift-and-invert mode, and basically find 416 // the largest eigenvalues of the inverse operator 417 // 418 mode = 3; 419 } 420 421 // The user-specified number of eigenvalues/vectors to compute 422 // 423 int nev = (int)nbrEigenvalues; 424 425 // Allocate space for ARPACK to store the residual 426 // 427 Scalar *resid = new Scalar[n]; 428 429 // Number of Lanczos vectors, must satisfy nev < ncv <= n 430 // Note that this indicates that nev != n, and we cannot compute 431 // all eigenvalues of a mtrix 432 // 433 int ncv = std::min(std::max(2*nev, 20), n); 434 435 // The working n x ncv matrix, also store the final eigenvectors (if computed) 436 // 437 Scalar *v = new Scalar[n*ncv]; 438 int ldv = n; 439 440 // Working space 441 // 442 Scalar *workd = new Scalar[3*n]; 443 int lworkl = ncv*ncv+8*ncv; // Must be at least this length 444 Scalar *workl = new Scalar[lworkl]; 445 446 int *iparam= new int[11]; 447 iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it 448 iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1))); 449 iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert 450 451 // Used during reverse communicate to notify where arrays start 452 // 453 int *ipntr = new int[11]; 454 455 // Error codes are returned in here, initial value of 0 indicates a random initial 456 // residual vector is used, any other values means resid contains the initial residual 457 // vector, possibly from a previous run 458 // 459 int info = 0; 460 461 Scalar scale = 1.0; 462 //if (!isBempty) 463 //{ 464 //Scalar scale = B.norm() / std::sqrt(n); 465 //scale = std::pow(2, std::floor(std::log(scale+1))); 466 ////M /= scale; 467 //for (size_t i=0; i<(size_t)B.outerSize(); i++) 468 // for (typename MatrixType::InnerIterator it(B, i); it; ++it) 469 // it.valueRef() /= scale; 470 //} 471 472 MatrixSolver OP; 473 if (mode == 1 || mode == 2) 474 { 475 if (!isBempty) 476 OP.compute(B); 477 } 478 else if (mode == 3) 479 { 480 if (sigma == 0.0) 481 { 482 OP.compute(A); 483 } 484 else 485 { 486 // Note: We will never enter here because sigma must be 0.0 487 // 488 if (isBempty) 489 { 490 MatrixType AminusSigmaB(A); 491 for (Index i=0; i<A.rows(); ++i) 492 AminusSigmaB.coeffRef(i,i) -= sigma; 493 494 OP.compute(AminusSigmaB); 495 } 496 else 497 { 498 MatrixType AminusSigmaB = A - sigma * B; 499 OP.compute(AminusSigmaB); 500 } 501 } 502 } 503 504 if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success) 505 std::cout << "Error factoring matrix" << std::endl; 506 507 do 508 { 509 internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, 510 &ncv, v, &ldv, iparam, ipntr, workd, workl, 511 &lworkl, &info); 512 513 if (ido == -1 || ido == 1) 514 { 515 Scalar *in = workd + ipntr[0] - 1; 516 Scalar *out = workd + ipntr[1] - 1; 517 518 if (ido == 1 && mode != 2) 519 { 520 Scalar *out2 = workd + ipntr[2] - 1; 521 if (isBempty || mode == 1) 522 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 523 else 524 Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 525 526 in = workd + ipntr[2] - 1; 527 } 528 529 if (mode == 1) 530 { 531 if (isBempty) 532 { 533 // OP = A 534 // 535 Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 536 } 537 else 538 { 539 // OP = L^{-1}AL^{-T} 540 // 541 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out); 542 } 543 } 544 else if (mode == 2) 545 { 546 if (ido == 1) 547 Matrix<Scalar, Dynamic, 1>::Map(in, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 548 549 // OP = B^{-1} A 550 // 551 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 552 } 553 else if (mode == 3) 554 { 555 // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I) 556 // The B * in is already computed and stored at in if ido == 1 557 // 558 if (ido == 1 || isBempty) 559 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 560 else 561 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n)); 562 } 563 } 564 else if (ido == 2) 565 { 566 Scalar *in = workd + ipntr[0] - 1; 567 Scalar *out = workd + ipntr[1] - 1; 568 569 if (isBempty || mode == 1) 570 Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 571 else 572 Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 573 } 574 } while (ido != 99); 575 576 if (info == 1) 577 m_info = NoConvergence; 578 else if (info == 3) 579 m_info = NumericalIssue; 580 else if (info < 0) 581 m_info = InvalidInput; 582 else if (info != 0) 583 eigen_assert(false && "Unknown ARPACK return value!"); 584 else 585 { 586 // Do we compute eigenvectors or not? 587 // 588 int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors; 589 590 // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK)) 591 // 592 char howmny[2] = "A"; 593 594 // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK) 595 // 596 int *select = new int[ncv]; 597 598 // Final eigenvalues 599 // 600 m_eivalues.resize(nev, 1); 601 602 internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv, 603 &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv, 604 v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info); 605 606 if (info == -14) 607 m_info = NoConvergence; 608 else if (info != 0) 609 m_info = InvalidInput; 610 else 611 { 612 if (rvec) 613 { 614 m_eivec.resize(A.rows(), nev); 615 for (int i=0; i<nev; i++) 616 for (int j=0; j<n; j++) 617 m_eivec(j,i) = v[i*n+j] / scale; 618 619 if (mode == 1 && !isBempty && BisSPD) 620 internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data()); 621 622 m_eigenvectorsOk = true; 623 } 624 625 m_nbrIterations = iparam[2]; 626 m_nbrConverged = iparam[4]; 627 628 m_info = Success; 629 } 630 631 delete[] select; 632 } 633 634 delete[] v; 635 delete[] iparam; 636 delete[] ipntr; 637 delete[] workd; 638 delete[] workl; 639 delete[] resid; 640 641 m_isInitialized = true; 642 643 return *this; 644 } 645 646 647 // Single precision 648 // 649 extern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which, 650 int *nev, float *tol, float *resid, int *ncv, 651 float *v, int *ldv, int *iparam, int *ipntr, 652 float *workd, float *workl, int *lworkl, 653 int *info); 654 655 extern "C" void sseupd_(int *rvec, char *All, int *select, float *d, 656 float *z, int *ldz, float *sigma, 657 char *bmat, int *n, char *which, int *nev, 658 float *tol, float *resid, int *ncv, float *v, 659 int *ldv, int *iparam, int *ipntr, float *workd, 660 float *workl, int *lworkl, int *ierr); 661 662 // Double precision 663 // 664 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which, 665 int *nev, double *tol, double *resid, int *ncv, 666 double *v, int *ldv, int *iparam, int *ipntr, 667 double *workd, double *workl, int *lworkl, 668 int *info); 669 670 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d, 671 double *z, int *ldz, double *sigma, 672 char *bmat, int *n, char *which, int *nev, 673 double *tol, double *resid, int *ncv, double *v, 674 int *ldv, int *iparam, int *ipntr, double *workd, 675 double *workl, int *lworkl, int *ierr); 676 677 678 namespace internal { 679 680 template<typename Scalar, typename RealScalar> struct arpack_wrapper 681 { 682 static inline void saupd(int *ido, char *bmat, int *n, char *which, 683 int *nev, RealScalar *tol, Scalar *resid, int *ncv, 684 Scalar *v, int *ldv, int *iparam, int *ipntr, 685 Scalar *workd, Scalar *workl, int *lworkl, int *info) 686 { 687 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 688 } 689 690 static inline void seupd(int *rvec, char *All, int *select, Scalar *d, 691 Scalar *z, int *ldz, RealScalar *sigma, 692 char *bmat, int *n, char *which, int *nev, 693 RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, 694 int *ldv, int *iparam, int *ipntr, Scalar *workd, 695 Scalar *workl, int *lworkl, int *ierr) 696 { 697 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 698 } 699 }; 700 701 template <> struct arpack_wrapper<float, float> 702 { 703 static inline void saupd(int *ido, char *bmat, int *n, char *which, 704 int *nev, float *tol, float *resid, int *ncv, 705 float *v, int *ldv, int *iparam, int *ipntr, 706 float *workd, float *workl, int *lworkl, int *info) 707 { 708 ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 709 } 710 711 static inline void seupd(int *rvec, char *All, int *select, float *d, 712 float *z, int *ldz, float *sigma, 713 char *bmat, int *n, char *which, int *nev, 714 float *tol, float *resid, int *ncv, float *v, 715 int *ldv, int *iparam, int *ipntr, float *workd, 716 float *workl, int *lworkl, int *ierr) 717 { 718 sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 719 workd, workl, lworkl, ierr); 720 } 721 }; 722 723 template <> struct arpack_wrapper<double, double> 724 { 725 static inline void saupd(int *ido, char *bmat, int *n, char *which, 726 int *nev, double *tol, double *resid, int *ncv, 727 double *v, int *ldv, int *iparam, int *ipntr, 728 double *workd, double *workl, int *lworkl, int *info) 729 { 730 dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 731 } 732 733 static inline void seupd(int *rvec, char *All, int *select, double *d, 734 double *z, int *ldz, double *sigma, 735 char *bmat, int *n, char *which, int *nev, 736 double *tol, double *resid, int *ncv, double *v, 737 int *ldv, int *iparam, int *ipntr, double *workd, 738 double *workl, int *lworkl, int *ierr) 739 { 740 dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 741 workd, workl, lworkl, ierr); 742 } 743 }; 744 745 746 template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> 747 struct OP 748 { 749 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out); 750 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs); 751 }; 752 753 template<typename MatrixSolver, typename MatrixType, typename Scalar> 754 struct OP<MatrixSolver, MatrixType, Scalar, true> 755 { 756 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 757 { 758 // OP = L^{-1} A L^{-T} (B = LL^T) 759 // 760 // First solve L^T out = in 761 // 762 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 763 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 764 765 // Then compute out = A out 766 // 767 Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n); 768 769 // Then solve L out = out 770 // 771 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 772 Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n)); 773 } 774 775 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 776 { 777 // Solve L^T out = in 778 // 779 Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k)); 780 Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k); 781 } 782 783 }; 784 785 template<typename MatrixSolver, typename MatrixType, typename Scalar> 786 struct OP<MatrixSolver, MatrixType, Scalar, false> 787 { 788 static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 789 { 790 eigen_assert(false && "Should never be in here..."); 791 } 792 793 static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 794 { 795 eigen_assert(false && "Should never be in here..."); 796 } 797 798 }; 799 800 } // end namespace internal 801 802 } // end namespace Eigen 803 804 #endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H 805 806
