17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 David Harmon <dharmon@gmail.com> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Eigen is free software; you can redistribute it and/or 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// modify it under the terms of the GNU Lesser General Public 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// License as published by the Free Software Foundation; either 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// version 3 of the License, or (at your option) any later version. 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Alternatively, you can redistribute it and/or 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// modify it under the terms of the GNU General Public License as 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// published by the Free Software Foundation; either version 2 of 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// the License, or (at your option) any later version. 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// GNU General Public License for more details. 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// You should have received a copy of the GNU Lesser General Public 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// License and a copy of the GNU General Public License along with 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Eigen. If not, see <http://www.gnu.org/licenses/>. 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_ARPACKGENERALIZEDSELFADJOINTEIGENSOLVER_H 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <Eigen/Dense> 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Scalar, typename RealScalar> struct arpack_wrapper; 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> struct OP; 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename MatrixSolver=SimplicialLLT<MatrixType>, bool BisSPD=false> 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass ArpackGeneralizedSelfAdjointEigenSolver 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //typedef typename MatrixSolver::MatrixType MatrixType; 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Scalar type for matrices of type \p MatrixType. */ 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Real scalar type for \p MatrixType. 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is just \c Scalar if #Scalar is real (e.g., \c float or 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \c Scalar), and the type of the real part of \c Scalar if #Scalar is 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * complex. 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename NumTraits<Scalar>::Real RealScalar; 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is a column vector with entries of type #RealScalar. 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The length of the vector is the size of \p nbrEigenvalues. 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default constructor. 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The default constructor is for cases in which the user intends to 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * perform decompositions via compute(). 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver() 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_eivec(), 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivalues(), 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false), 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eigenvectorsOk(false), 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrConverged(0), 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrIterations(0) 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { } 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Constructor; computes generalized eigenvalues of given matrix with respect to another matrix. 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * computed. By default, the upper triangular part is used, but can be changed 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * through the template parameter. 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Self-adjoint matrix for the generalized eigenvalue problem. 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Must be less than the size of the input matrix, or an error is returned. 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * respective meanings to find the largest magnitude , smallest magnitude, 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * value can contain floating point value in string form, in which case the 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues closest to this value will be found. 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * means machine precision. 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This constructor calls compute(const MatrixType&, const MatrixType&, Index, string, int, RealScalar) 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * to compute the eigenvalues of the matrix \p A with respect to \p B. The eigenvectors are computed if 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \p options equals #ComputeEigenvectors. 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, const MatrixType& B, 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nbrEigenvalues, std::string eigs_sigma="LM", 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int options=ComputeEigenvectors, RealScalar tol=0.0) 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_eivec(), 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivalues(), 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false), 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eigenvectorsOk(false), 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrConverged(0), 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrIterations(0) 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Constructor; computes eigenvalues of given matrix. 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Self-adjoint matrix whose eigenvalues / eigenvectors will 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * computed. By default, the upper triangular part is used, but can be changed 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * through the template parameter. 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Must be less than the size of the input matrix, or an error is returned. 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * respective meanings to find the largest magnitude , smallest magnitude, 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * value can contain floating point value in string form, in which case the 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues closest to this value will be found. 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * means machine precision. 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This constructor calls compute(const MatrixType&, Index, string, int, RealScalar) 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * to compute the eigenvalues of the matrix \p A. The eigenvectors are computed if 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \p options equals #ComputeEigenvectors. 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver(const MatrixType& A, 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nbrEigenvalues, std::string eigs_sigma="LM", 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int options=ComputeEigenvectors, RealScalar tol=0.0) 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_eivec(), 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivalues(), 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false), 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eigenvectorsOk(false), 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrConverged(0), 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrIterations(0) 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(A, nbrEigenvalues, eigs_sigma, options, tol); 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes generalized eigenvalues / eigenvectors of given matrix using the external ARPACK library. 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Selfadjoint matrix for generalized eigenvalues. 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Must be less than the size of the input matrix, or an error is returned. 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * respective meanings to find the largest magnitude , smallest magnitude, 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * value can contain floating point value in string form, in which case the 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues closest to this value will be found. 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * means machine precision. 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns Reference to \c *this 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This function computes the generalized eigenvalues of \p A with respect to \p B using ARPACK. The eigenvalues() 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * then the eigenvectors are also computed and can be retrieved by 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * calling eigenvectors(). 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, const MatrixType& B, 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nbrEigenvalues, std::string eigs_sigma="LM", 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int options=ComputeEigenvectors, RealScalar tol=0.0); 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes eigenvalues / eigenvectors of given matrix using the external ARPACK library. 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Selfadjoint matrix whose eigendecomposition is to be computed. 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] nbrEigenvalues The number of eigenvalues / eigenvectors to compute. 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Must be less than the size of the input matrix, or an error is returned. 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] eigs_sigma String containing either "LM", "SM", "LA", or "SA", with 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * respective meanings to find the largest magnitude , smallest magnitude, 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * largest algebraic, or smallest algebraic eigenvalues. Alternatively, this 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * value can contain floating point value in string form, in which case the 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues closest to this value will be found. 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] tol What tolerance to find the eigenvalues to. Default is 0, which 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * means machine precision. 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns Reference to \c *this 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This function computes the eigenvalues of \p A using ARPACK. The eigenvalues() 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * then the eigenvectors are also computed and can be retrieved by 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * calling eigenvectors(). 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver& compute(const MatrixType& A, 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nbrEigenvalues, std::string eigs_sigma="LM", 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int options=ComputeEigenvectors, RealScalar tol=0.0); 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the eigenvectors of given matrix. 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the matrix whose columns are the eigenvectors. 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre The eigenvectors have been computed before. 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvectors are normalized to have (Euclidean) norm equal to one. If 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * this object was used to solve the eigenproblem for the selfadjoint 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrix \f$ A \f$, then the matrix returned by this function is the 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrix \f$ V \f$ in the eigendecomposition \f$ A V = D V \f$. 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * For the generalized eigenproblem, the matrix returned is the solution \f$ A V = D B V \f$ 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa eigenvalues() 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Matrix<Scalar, Dynamic, Dynamic>& eigenvectors() const 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_eivec; 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the eigenvalues of given matrix. 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the column vector containing the eigenvalues. 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre The eigenvalues have been computed before. 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The eigenvalues are repeated according to their algebraic multiplicity, 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * so there are as many eigenvalues as rows in the matrix. The eigenvalues 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * are sorted in increasing order. 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa eigenvectors(), MatrixBase::eigenvalues() 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Matrix<Scalar, Dynamic, 1>& eigenvalues() const 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_eivalues; 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes the positive-definite square root of the matrix. 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns the positive-definite square root of the matrix 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre The eigenvalues and eigenvectors of a positive-definite matrix 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * have been computed before. 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The square root of a positive-definite matrix \f$ A \f$ is the 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * positive-definite matrix whose square equals \f$ A \f$. This function 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$. 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa operatorInverseSqrt(), 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \ref MatrixFunctions_Module "MatrixFunctions Module" 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Dynamic> operatorSqrt() const 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes the inverse square root of the matrix. 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns the inverse positive-definite square root of the matrix 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre The eigenvalues and eigenvectors of a positive-definite matrix 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * have been computed before. 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * cheaper than first computing the square root with operatorSqrt() and 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * then its inverse with MatrixBase::inverse(). 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa operatorSqrt(), MatrixBase::inverse(), 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \ref MatrixFunctions_Module "MatrixFunctions Module" 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Dynamic> operatorInverseSqrt() const 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Reports whether previous computation was successful. 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo info() const 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "ArpackGeneralizedSelfAdjointEigenSolver is not initialized."); 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_info; 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez size_t getNbrConvergedEigenValues() const 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { return m_nbrConverged; } 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez size_t getNbrIterations() const 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { return m_nbrIterations; } 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprotected: 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Dynamic> m_eivec; 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1> m_eivalues; 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo m_info; 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_isInitialized; 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_eigenvectorsOk; 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez size_t m_nbrConverged; 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez size_t m_nbrIterations; 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename MatrixSolver, bool BisSPD> 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez::compute(const MatrixType& A, Index nbrEigenvalues, 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::string eigs_sigma, int options, RealScalar tol) 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType B(0,0); 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(A, B, nbrEigenvalues, eigs_sigma, options, tol); 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename MatrixSolver, bool BisSPD> 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD>& 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ArpackGeneralizedSelfAdjointEigenSolver<MatrixType, MatrixSolver, BisSPD> 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez::compute(const MatrixType& A, const MatrixType& B, Index nbrEigenvalues, 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::string eigs_sigma, int options, RealScalar tol) 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(A.cols() == A.rows()); 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(B.cols() == B.rows()); 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(B.rows() == 0 || A.cols() == B.rows()); 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert((options &~ (EigVecMask | GenEigMask)) == 0 3557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && (options & EigVecMask) != EigVecMask 3567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && "invalid option parameter"); 3577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool isBempty = (B.rows() == 0) || (B.cols() == 0); 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // For clarity, all parameters match their ARPACK name 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Always 0 on the first call 3637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int ido = 0; 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int n = (int)A.cols(); 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // User options: "LA", "SA", "SM", "LM", "BE" 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char whch[3] = "LM"; 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Specifies the shift if iparam[6] = { 3, 4, 5 }, not used if iparam[6] = { 1, 2 } 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar sigma = 0.0; 3757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (eigs_sigma.length() >= 2 && isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigs_sigma[0] = toupper(eigs_sigma[0]); 3797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigs_sigma[1] = toupper(eigs_sigma[1]); 3807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // In the following special case we're going to invert the problem, since solving 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // for larger magnitude is much much faster 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // i.e., if 'SM' is specified, we're going to really use 'LM', the default 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (eigs_sigma.substr(0,2) != "SM") 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez whch[0] = eigs_sigma[0]; 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez whch[1] = eigs_sigma[1]; 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false && "Specifying clustered eigenvalues is not yet supported!"); 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If it's not scalar values, then the user may be explicitly 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // specifying the sigma value to cluster the evs around 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez sigma = atof(eigs_sigma.c_str()); 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If atof fails, it returns 0.0, which is a fine default 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // "I" means normal eigenvalue problem, "G" means generalized 4057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char bmat[2] = "I"; 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1])) || (!isBempty && !BisSPD)) 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bmat[0] = 'G'; 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Now we determine the mode to use 4117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int mode = (bmat[0] == 'G') + 1; 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (eigs_sigma.substr(0,2) == "SM" || !(isalpha(eigs_sigma[0]) && isalpha(eigs_sigma[1]))) 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // We're going to use shift-and-invert mode, and basically find 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // the largest eigenvalues of the inverse operator 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mode = 3; 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // The user-specified number of eigenvalues/vectors to compute 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int nev = (int)nbrEigenvalues; 4247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Allocate space for ARPACK to store the residual 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *resid = new Scalar[n]; 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Number of Lanczos vectors, must satisfy nev < ncv <= n 4307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Note that this indicates that nev != n, and we cannot compute 4317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // all eigenvalues of a mtrix 4327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int ncv = std::min(std::max(2*nev, 20), n); 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // The working n x ncv matrix, also store the final eigenvectors (if computed) 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *v = new Scalar[n*ncv]; 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int ldv = n; 4397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Working space 4417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *workd = new Scalar[3*n]; 4437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int lworkl = ncv*ncv+8*ncv; // Must be at least this length 4447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *workl = new Scalar[lworkl]; 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *iparam= new int[11]; 4477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez iparam[0] = 1; // 1 means we let ARPACK perform the shifts, 0 means we'd have to do it 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez iparam[2] = std::max(300, (int)std::ceil(2*n/std::max(ncv,1))); 4497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez iparam[6] = mode; // The mode, 1 is standard ev problem, 2 for generalized ev, 3 for shift-and-invert 4507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Used during reverse communicate to notify where arrays start 4527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ipntr = new int[11]; 4547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Error codes are returned in here, initial value of 0 indicates a random initial 4567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // residual vector is used, any other values means resid contains the initial residual 4577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // vector, possibly from a previous run 4587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info = 0; 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar scale = 1.0; 4627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //if (!isBempty) 4637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //{ 4647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Scalar scale = B.norm() / std::sqrt(n); 4657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //scale = std::pow(2, std::floor(std::log(scale+1))); 4667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ////M /= scale; 4677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //for (size_t i=0; i<(size_t)B.outerSize(); i++) 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // for (typename MatrixType::InnerIterator it(B, i); it; ++it) 4697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // it.valueRef() /= scale; 4707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //} 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixSolver OP; 4737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (mode == 1 || mode == 2) 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (!isBempty) 4767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OP.compute(B); 4777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (mode == 3) 4797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (sigma == 0.0) 4817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OP.compute(A); 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Note: We will never enter here because sigma must be 0.0 4877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (isBempty) 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType AminusSigmaB(A); 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i=0; i<A.rows(); ++i) 4927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez AminusSigmaB.coeffRef(i,i) -= sigma; 4937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OP.compute(AminusSigmaB); 4957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType AminusSigmaB = A - sigma * B; 4997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OP.compute(AminusSigmaB); 5007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (!(mode == 1 && isBempty) && !(mode == 2 && isBempty) && OP.info() != Success) 5057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cout << "Error factoring matrix" << std::endl; 5067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez do 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::arpack_wrapper<Scalar, RealScalar>::saupd(&ido, bmat, &n, whch, &nev, &tol, resid, 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez &ncv, v, &ldv, iparam, ipntr, workd, workl, 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez &lworkl, &info); 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (ido == -1 || ido == 1) 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *in = workd + ipntr[0] - 1; 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *out = workd + ipntr[1] - 1; 5177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (ido == 1 && mode != 2) 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *out2 = workd + ipntr[2] - 1; 5217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (isBempty || mode == 1) 5227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out2, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out2, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 5257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez in = workd + ipntr[2] - 1; 5277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (mode == 1) 5307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (isBempty) 5327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // OP = A 5347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 5367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // OP = L^{-1}AL^{-T} 5407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::applyOP(OP, A, n, in, out); 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (mode == 2) 5457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (ido == 1) 5477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(in, n) = A * Matrix<Scalar, Dynamic, 1>::Map(in, n); 5487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // OP = B^{-1} A 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (mode == 3) 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // OP = (A-\sigmaB)B (\sigma could be 0, and B could be I) 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // The B * in is already computed and stored at in if ido == 1 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (ido == 1 || isBempty) 5597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 5607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.solve(B * Matrix<Scalar, Dynamic, 1>::Map(in, n)); 5627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (ido == 2) 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *in = workd + ipntr[0] - 1; 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *out = workd + ipntr[1] - 1; 5687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (isBempty || mode == 1) 5707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = Matrix<Scalar, Dynamic, 1>::Map(in, n); 5717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = B * Matrix<Scalar, Dynamic, 1>::Map(in, n); 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } while (ido != 99); 5757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (info == 1) 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NoConvergence; 5787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (info == 3) 5797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (info < 0) 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = InvalidInput; 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (info != 0) 5837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false && "Unknown ARPACK return value!"); 5847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Do we compute eigenvectors or not? 5877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int rvec = (options & ComputeEigenvectors) == ComputeEigenvectors; 5897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // "A" means "All", use "S" to choose specific eigenvalues (not yet supported in ARPACK)) 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char howmny[2] = "A"; 5937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // if howmny == "S", specifies the eigenvalues to compute (not implemented in ARPACK) 5957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 5967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *select = new int[ncv]; 5977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Final eigenvalues 5997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 6007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivalues.resize(nev, 1); 6017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::arpack_wrapper<Scalar, RealScalar>::seupd(&rvec, howmny, select, m_eivalues.data(), v, &ldv, 6037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez &sigma, bmat, &n, whch, &nev, &tol, resid, &ncv, 6047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez v, &ldv, iparam, ipntr, workd, workl, &lworkl, &info); 6057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (info == -14) 6077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NoConvergence; 6087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (info != 0) 6097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = InvalidInput; 6107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 6117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (rvec) 6137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivec.resize(A.rows(), nev); 6157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (int i=0; i<nev; i++) 6167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (int j=0; j<n; j++) 6177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eivec(j,i) = v[i*n+j] / scale; 6187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (mode == 1 && !isBempty && BisSPD) 6207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::OP<MatrixSolver, MatrixType, Scalar, BisSPD>::project(OP, n, nev, m_eivec.data()); 6217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_eigenvectorsOk = true; 6237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrIterations = iparam[2]; 6267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_nbrConverged = iparam[4]; 6277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = Success; 6297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete select; 6327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete v; 6357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete iparam; 6367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete ipntr; 6377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete workd; 6387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete workl; 6397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez delete resid; 6407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized = true; 6427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 6447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 6457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Single precision 6487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 6497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezextern "C" void ssaupd_(int *ido, char *bmat, int *n, char *which, 6507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *nev, float *tol, float *resid, int *ncv, 6517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *v, int *ldv, int *iparam, int *ipntr, 6527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *workd, float *workl, int *lworkl, 6537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *info); 6547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezextern "C" void sseupd_(int *rvec, char *All, int *select, float *d, 6567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *z, int *ldz, float *sigma, 6577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char *bmat, int *n, char *which, int *nev, 6587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *tol, float *resid, int *ncv, float *v, 6597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ldv, int *iparam, int *ipntr, float *workd, 6607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *workl, int *lworkl, int *ierr); 6617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Double precision 6637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 6647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezextern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which, 6657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *nev, double *tol, double *resid, int *ncv, 6667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *v, int *ldv, int *iparam, int *ipntr, 6677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *workd, double *workl, int *lworkl, 6687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *info); 6697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezextern "C" void dseupd_(int *rvec, char *All, int *select, double *d, 6717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *z, int *ldz, double *sigma, 6727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char *bmat, int *n, char *which, int *nev, 6737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *tol, double *resid, int *ncv, double *v, 6747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ldv, int *iparam, int *ipntr, double *workd, 6757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *workl, int *lworkl, int *ierr); 6767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 6797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar, typename RealScalar> struct arpack_wrapper 6817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void saupd(int *ido, char *bmat, int *n, char *which, 6837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *nev, RealScalar *tol, Scalar *resid, int *ncv, 6847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *v, int *ldv, int *iparam, int *ipntr, 6857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *workd, Scalar *workl, int *lworkl, int *info) 6867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 6887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void seupd(int *rvec, char *All, int *select, Scalar *d, 6917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *z, int *ldz, RealScalar *sigma, 6927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char *bmat, int *n, char *which, int *nev, 6937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar *tol, Scalar *resid, int *ncv, Scalar *v, 6947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ldv, int *iparam, int *ipntr, Scalar *workd, 6957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar *workl, int *lworkl, int *ierr) 6967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) 6987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <> struct arpack_wrapper<float, float> 7027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void saupd(int *ido, char *bmat, int *n, char *which, 7047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *nev, float *tol, float *resid, int *ncv, 7057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *v, int *ldv, int *iparam, int *ipntr, 7067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *workd, float *workl, int *lworkl, int *info) 7077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ssaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 7097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void seupd(int *rvec, char *All, int *select, float *d, 7127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *z, int *ldz, float *sigma, 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char *bmat, int *n, char *which, int *nev, 7147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *tol, float *resid, int *ncv, float *v, 7157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ldv, int *iparam, int *ipntr, float *workd, 7167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez float *workl, int *lworkl, int *ierr) 7177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez sseupd_(rvec, All, select, d, z, ldz, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 7197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez workd, workl, lworkl, ierr); 7207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <> struct arpack_wrapper<double, double> 7247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void saupd(int *ido, char *bmat, int *n, char *which, 7267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *nev, double *tol, double *resid, int *ncv, 7277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *v, int *ldv, int *iparam, int *ipntr, 7287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *workd, double *workl, int *lworkl, int *info) 7297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dsaupd_(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info); 7317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void seupd(int *rvec, char *All, int *select, double *d, 7347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *z, int *ldz, double *sigma, 7357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez char *bmat, int *n, char *which, int *nev, 7367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *tol, double *resid, int *ncv, double *v, 7377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int *ldv, int *iparam, int *ipntr, double *workd, 7387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double *workl, int *lworkl, int *ierr) 7397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dseupd_(rvec, All, select, d, v, ldv, sigma, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, 7417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez workd, workl, lworkl, ierr); 7427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixSolver, typename MatrixType, typename Scalar, bool BisSPD> 7477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct OP 7487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out); 7507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs); 7517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixSolver, typename MatrixType, typename Scalar> 7547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct OP<MatrixSolver, MatrixType, Scalar, true> 7557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 7577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // OP = L^{-1} A L^{-T} (B = LL^T) 7597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 7607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // First solve L^T out = in 7617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 7627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixU().solve(Matrix<Scalar, Dynamic, 1>::Map(in, n)); 7637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationPinv() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 7647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Then compute out = A out 7667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 7677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = A * Matrix<Scalar, Dynamic, 1>::Map(out, n); 7687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Then solve L out = out 7707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 7717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.permutationP() * Matrix<Scalar, Dynamic, 1>::Map(out, n); 7727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, 1>::Map(out, n) = OP.matrixL().solve(Matrix<Scalar, Dynamic, 1>::Map(out, n)); 7737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 7767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Solve L^T out = in 7787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // 7797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.matrixU().solve(Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k)); 7807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k) = OP.permutationPinv() * Matrix<Scalar, Dynamic, Dynamic>::Map(vecs, n, k); 7817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixSolver, typename MatrixType, typename Scalar> 7867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct OP<MatrixSolver, MatrixType, Scalar, false> 7877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void applyOP(MatrixSolver &OP, const MatrixType &A, int n, Scalar *in, Scalar *out) 7897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false && "Should never be in here..."); 7917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static inline void project(MatrixSolver &OP, int n, int k, Scalar *vecs) 7947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false && "Should never be in here..."); 7967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 8017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 8037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_ARPACKSELFADJOINTEIGENSOLVER_H 8057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 806