1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{{\libtitle eigenCore} sub-library} 18 19Most of the standard eigenvalue algorithms exploit projection processes in order to extract approximate eigenvectors from a given subspace. The basic idea of a projection method is to extract an approximate eigenvector from a specified low-dimensional subspace. If these approximation can be extracted, a small matrix eigenvalue problem is obtained. The numerical solution of the small $m\times m$ eigenvalue problem will be treated by the {\lib eigenCore} library. 20The {\lib eigenCore} library is organized in multi-level, each of which utilizes subroutines provided by the lower one. 21 \begin{itemize} 22\item {\lib utils} : the lowest level provides various basic operations 23\item {\lib decomposition} : the second level implements several decomposition algorithms 24\item {\lib eigenSolver} : the highest level solves different eigen problem 25\end{itemize} 26 27\subsection{{\libtitle utils} part} 28 29This part contains several classes which work on different basic operations. 30 31\subsubsection{The {\classtitle VectorEigenDense} class} 32 33Class {\class VectorEigenDense} is derived directly from {\class Vector}. Besides taking advantages of all funtionalities of {\class Vector}, this class supplements some particular functions to calculate Householder transformation. 34\vspace{.1cm} 35\begin{lstlisting} 36template<typename K> 37class VectorEigenDense : public Vector<K> 38{ 39public: 40 typedef typename VectorEigenDense::type_t Scalar; 41 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 42 typedef typename VectorEigenDense::it_vk it_vk; 43 typedef typename VectorEigenDense::cit_vk cit_vk; 44 45 VectorEigenDense() : Vector<K>(), acType_(_col) {} 46 VectorEigenDense(const Dimen l) : Vector<K>(l, K()), acType_(_col) {} 47 ... 48 49\end{lstlisting} 50\vspace{.1cm} 51This class offers some additional functions to calculate House holder values. 52\vspace{.1cm} 53\begin{lstlisting}[]{} 54void makeHouseHolderInPlace(K& tau, Real& beta); 55void makeHouseHolderInPlace(K& tau, Real& beta, Number tail); 56void makeHouseHolder(VectorEigenDense<K>& essential, K& tau, Real& beta) const; 57\end{lstlisting} 58\vspace{.2cm} 59Note that the in the future, we may improve calculation performance of eigensolver for {\em dense} matrix by changing this class {\class VectorEigenDense}. By using the same block of allocated data, we can gain more efficiencies. However, for the moment, it's sufficient to have this class inherited from class {\class Vector} 60\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}} 61 62\subsubsection{The {\classtitle MatrixEigenDense} class} 63 64Class {\class MatrixEigenDense} inherits from the fundamental class {\class Matrix}. Basically, this class reuses all functionalities of its father in addition to some methods to support computing House holder, Cholesky composition and Jacobian rotation. 65\vspace{.1cm} 66\begin{lstlisting} 67template<typename K> 68class MatrixEigenDense : public Matrix<K> 69{ 70public: 71 typedef typename MatrixEigenDense<K>::type_t Scalar; 72 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 73 typedef typename std::vector<K>::iterator it_vk; 74 typedef typename std::vector<K>::const_iterator cit_vk; 75 76 MatrixEigenDense() : Matrix<K>(), col_(1) {} 77 MatrixEigenDense(const Dimen r) : Matrix<K>(r,r), col_(r) {} 78 ... 79 80\end{lstlisting} 81\vspace{.1cm} 82Householder computation can be invoked with 83\vspace{.1cm} 84\begin{lstlisting}[]{} 85void applyHouseholderOnTheLeft(const VectorEigenDense<K>& essential, const K& tau); 86void applyHouseholderOnTheRight(const VectorEigenDense<K>& essential, const K& tau); 87\end{lstlisting} 88\vspace{.1cm} 89Cholesky decomposition can be calculated by 90\vspace{.1cm} 91\begin{lstlisting}[]{} 92MatrixEigenDense<K> cholesky() const; 93\end{lstlisting} 94\vspace{.1cm} 95or one can call Jacobian rotation with 96\vspace{.1cm} 97\begin{lstlisting}[]{} 98template<typename OtherScalar> 99inline void applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j); 100template<typename OtherScalar> 101inline void applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j); 102\end{lstlisting} 103\vspace{.1cm} 104Similar to class {\class VectorEigenDense}, in the future, this class may be improved to make {\lib eigenCore} more efficient. 105\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}} 106 107\subsubsection{The {\classtitle JacobiRotation} class} 108 109Class {\class JacobiRotation} represents a {\em Jacobian} or {\em Givens} rotation, 2D rotation in a plane defined by its cosine and sine. 110\vspace{.1cm} 111\begin{lstlisting} 112template<typename Scalar> 113class JacobiRotation 114{ 115 public: 116 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 117 118 /** Default constructor without any initialization. */ 119 JacobiRotation() {} 120 121 /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */ 122 JacobiRotation(const Scalar& c, const Scalar& s) : c_(c), s_(s) {} 123 124 Scalar& c() { return c_; } 125 Scalar c() const { return c_; } 126 Scalar& s() { return s_; } 127 Scalar s() const { return s_; } 128 ... 129\end{lstlisting} 130\vspace{.1cm} 131The class provide two basic operations: making a Jacobi rotation and applying this rotation on both the right and left side of a selfadjoint $2\times2$ matrix; making a Givens rotation and applying this rotation on the left side of a vector. 132\vspace{.1cm} 133\begin{lstlisting}[deletekeywords={[3] x, y, z}] 134template<typename Scalar> 135bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z); 136template<typename Scalar> 137void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z); 138\end{lstlisting} 139\vspace{.1cm} 140 141\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={Traits.hpp}} 142 143\subsubsection{The {\classtitle NumTraits} class} 144 145{\class NumTraits} greatly facilitate the management of the sort of extra parameters that come up during the implementation of eigensolver algorithm. 146\begin{lstlisting} 147template<typename K> 148struct NumTraits {}; 149 150template<> 151struct NumTraits<Real> { 152 enum { 153 IsComplex = 0 154 }; 155 static inline bool isReal() { return true;} 156 static inline bool isComplex() { return false;} 157 static inline Real zero() { return Real(0.0);} 158 static inline Real one() { return Real(1.0);} 159 static inline Real imag(const Real& v) { return v;} 160 static inline Real real(const Real& v) { return v;} 161 static inline Real value(const Real& v) { return v;} 162 static inline Real abs2(const Real& v) { return v*v;} 163 static inline Real epsilon() { return std::numeric_limits<Real>::epsilon(); } 164 typedef Real RealScalar; 165 typedef Complex ComplexScalar; 166}; 167 168template<> 169struct NumTraits<Complex> { 170 enum { 171 IsComplex = 1 172 }; 173 static inline bool isReal() {return false;} 174 static inline bool isComplex() { return true;} 175 static inline Complex zero() { return Complex(0.0,0.0);} 176 static inline Complex one() { return Complex(1.0,0.0);} 177 static inline Real imag(const Complex& v) { return v.imag();} 178 static inline Real real(const Complex& v) { return v.real();} 179 static inline Complex value(const Complex& v) { return v;} 180 static inline Real abs2(const Complex& v) { return v.real()*v.real() + v.imag()*v.imag();} 181 typedef Real RealScalar; 182 typedef Complex ComplexScalar; 183}; 184\end{lstlisting} 185\vspace{.1cm} 186\displayInfos{library=eigenSolvers, test={test\_EigenSolverUtils.cpp}, header dep={utils.h}} 187 188\subsection{{\libtitle decomposition} part} 189 190All classes of this section utilize functions provided from the lower level - {\class utils} - to implement various decomposition algorithms, which play an important role in solving different eigen problems. 191 192\subsubsection{The {\classtitle Tridiagonalization} class} 193 194Class {\class Tridiagonalization} performs a tridiagonal decomposition of a self-adjoint matrix $ A$ such that: 195 $A = Q T Q^* $ where $Q$ is unitary and $ T $ a real symmetric tridiagonal matrix. A tridiagonal matrix is a matrix which has non-zero elements only on the main diagonal and the first diagonal below and above it. This class is used in SelfAdjointEigenSolver to compute the eigenvalues and eigenvectors of a selfadjoint matrix. 196\vspace{.1cm} 197\begin{lstlisting}[deletekeywords={[3] size}] 198template<typename _MatrixType> 199class Tridiagonalization 200{ 201 public: 202 203 /** \brief Synonym for the template parameter \p _MatrixType. */ 204 typedef _MatrixType MatrixType; 205 206 typedef typename MatrixType::type_t Scalar; 207 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 208 typedef VectorEigenDense<Scalar> CoeffVectorType; 209 typedef VectorEigenDense<RealScalar> DiagonalType; 210 ... 211 Tridiagonalization(Dimen size); 212 Tridiagonalization(const MatrixType& matrix); 213 ... 214\end{lstlisting} 215\vspace{.1cm} 216This class comes with default constructor with size of the matrix whose tridiagonal decomposition will be computed. In fact, this parameter is just a hint. After being created, object of this class can be used to calculate any square {\em dense} matrix with 217\vspace{.1cm} 218\begin{lstlisting}[]{} 219Tridiagonalization& compute(const MatrixType& matrix); 220\end{lstlisting} 221\vspace{.1cm} 222The unitary matrix $Q$ can only be invoked after calling the {\itshape compute}: 223\vspace{.1cm} 224\begin{lstlisting}[]{} 225MatrixType matrixQ(); 226\end{lstlisting} 227\vspace{.1cm} 228"Householder coefficients" can be retrieved by 229\vspace{.1cm} 230\begin{lstlisting}[]{} 231inline CoeffVectorType householderCoefficients() const; 232\end{lstlisting} 233\vspace{.1cm} 234We don't need the full-form of matrix $T$ in computation of eigenvalues; however its real diagonal and its real sub-diagonal are essential in calculating eigenvalues of the matrix. These two values can be retrieved with: 235\vspace{.1cm} 236\begin{lstlisting}[]{} 237DiagonalReturnType diagonal() const; 238SubDiagonalReturnType subDiagonal() const; 239\end{lstlisting} 240\vspace{.1cm} 241For example, the following $5\times 5$ real symmetric matrix 242\[ 243A=\left[ 244\begin{array}{cccccc} 2451.36 & -0.816 & 0.521 & 1.43 & -0.144\\ 246-0.816 & -0.659 & 0.794 & -0.173 & -0.406 \\ 2470.521 & 0.794 & -0.541 & 0.461& 0.179 \\ 2481.43& -0.173 & 0.461& -1.43 & 0.822 \\ 249-0.144 & -0.406 & 0.179 & 0.822& -1.37 250\end{array} 251\right] 252\] 253after being tridiagonalized, has diagonal and sub-diagonal vectors\\ 254\mbox{}\hspace{10 mm} diagonal = ( 1.36 -0.659 -0.541 -1.43 -1.37 ) \\ 255\mbox{}\hspace{10 mm} sub-diagonal = ( -0.816 0.794 0.461 0.822 ) 256 257\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Traits.hpp, VectorEigenDense.hpp, HouseHolderSequence.hpp}} 258 259\subsubsection{The {\classtitle HessenbergDecomposition} class } 260 261Class {\class HessenbergDecomposition} performs an Hessenberg decomposition of a matrix $ A $. In the real case, the Hessenberg decomposition consists of an orthogonal matrix $ Q $ and a Hessenberg matrix $ H $ such that $ A = Q H Q^T $. An orthogonal matrix is a matrix whose inverse equals its transpose ($ Q^{-1} = Q^T $). A Hessenberg matrix has zeros below the subdiagonal, so it is almost upper triangular. The Hessenberg decomposition of a complex matrix is $ A = Q H Q^* $ with $ Q $ unitary (that is, $ Q^{-1} = Q^* $). \\ 262\vspace{.1cm} 263\begin{lstlisting}[deletekeywords={[3] size}] 264template<typename _MatrixType> 265class HessenbergDecomposition 266{ 267 public: 268 269 /** \brief Synonym for the template parameter \p _MatrixType. */ 270 typedef _MatrixType MatrixType; 271 272 /** \brief Scalar type for matrices of type #MatrixType. */ 273 typedef typename MatrixType::type_t Scalar; 274 275 /** \brief Type for vector of Householder coefficients. 276 * 277 * This is column vector with entries of type #Scalar. The length of the 278 * vector is one less than the size of #MatrixType, 279 */ 280 typedef VectorEigenDense<Scalar> CoeffVectorType; 281 282 /** \brief Return type of matrixQ() */ 283 typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType; 284 285 /** \brief Default constructor; the decomposition will be computed later. 286 * 287 * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed. 288 * 289 * The default constructor is useful in cases in which the user intends to 290 * perform decompositions via compute(). The \p size parameter is only 291 * used as a hint. It is not an error to give a wrong \p size, but it may 292 * impair performance. 293 * 294 * \sa compute() for an example. 295 */ 296 HessenbergDecomposition(Number size) 297 : matrix_(size,size), 298 isInitialized_(false) 299 { 300 if(size>1) 301 hCoeffs_.resize(size-1); 302 } 303 ... 304\end{lstlisting} 305\vspace{.1cm} 306We can use the constructor which computes the Hessenberg decomposition at construction time. 307\vspace{.1cm} 308\begin{lstlisting}[]{} 309HessenbergDecomposition(const MatrixType&); 310\end{lstlisting} 311\vspace{.1cm} 312Once the decomposition is computed, we can use the {\cmd matrixH()} and {\cmd matrixQ()} functions to construct the matrices $H$ and $Q$ in the decomposition. 313\vspace{.1cm} 314\begin{lstlisting}[]{} 315MatrixType matrixQ(); 316MatrixType matrixH() const; 317\end{lstlisting} 318\vspace{.1cm} 319Alternatively, calling the function {\cmd compute()} to compute the Hessenberg decomposition of a given matrix. then we can use these two functions above 320\vspace{.1cm} 321\begin{lstlisting}[]{} 322HessenbergDecomposition& compute(const MatrixType& matrix); 323\end{lstlisting} 324\vspace{.1cm} 325For instance, the following $5\times 5$ real symmetric matrix 326\[ 327A=\left[ 328\begin{array}{cccccc} 3291.36 & -0.816 & 0.521 & 1.43 & -0.144\\ 330-0.816 & -0.659 & 0.794 & -0.173 & -0.406 \\ 3310.521 & 0.794 & -0.541 & 0.461& 0.179 \\ 3321.43& -0.173 & 0.461& -1.43 & 0.822 \\ 333-0.144 & -0.406 & 0.179 & 0.822& -1.37 334\end{array} 335\right] 336\] 337has Hessenberg decomposition 338\footnotesize 339\[ 340H=\left[ 341\begin{array}{cccccc} 3421.36 & 1.7328972849 & 0 & 0 & 0 \\ 3431.7328972849& -1.1933420269& -0.96533126903& -9.1846580094e-17& 2.021583842e-16 \\ 3440 & -0.96533126903& -1.2812029745& 0.21380585197& 1.9428902931e-16 \\ 3450 & 0 & 0.21380585197 & -1.6885671008& 0.34719629211 \\ 3460 & 0 & 0& 0.34719629211& 0.16311210222 347\end{array} 348\right] 349\] 350\normalsize 351and unitary matrix 352\footnotesize 353\[ 354Q=\left[ 355\begin{array}{cccccc} 3561 & 0 & 0 & 0 & 0\\ 3570& -0.47088769029& 0.12629876626& -0.67227759101& -0.55709626223 \\ 3580& 0.30065255716& -0.19453390601& 0.43672469898 & -0.82524913607 \\ 3590& 0.82520759451& 0.045097442878& -0.56297063107& -0.0079192904815 \\ 3600 & -0.083097827699& -0.97168482632& -0.20094388887& 0.092438643767 361\end{array} 362\right] 363\] 364\normalsize 365 366\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Traits.hpp, VectorEigenDense.hpp}} 367 368\subsubsection{The {\classtitle RealSchur} class } 369 370Class {\class RealSchur} performs a real Schur decomposition of a square matrix. Given a real square matrix $A$, this class computes the real Schur decomposition: $ A = U T U^T $ where $U$ is a real orthogonal matrix and $T$ is a real quasi-triangular matrix. An orthogonal matrix is a matrix whose inverse is equal to its transpose, $ U^{-1} = U^T $. A quasi-triangular matrix is a block-triangular matrix whose diagonal consists of 1-by-1 blocks and 2-by-2 blocks with complex eigenvalues. The eigenvalues of the blocks on the diagonal of $T$ are the same as the eigenvalues of the matrix $A$, and thus the real Schur decomposition is used in {\class RealEigenSolver} to compute the eigendecomposition of a matrix. 371\vspace{.1cm} 372\begin{lstlisting}[deletekeywords={[3] size}] 373template<typename _MatrixType> 374class RealSchur 375{ 376 public: 377 typedef _MatrixType MatrixType; 378 typedef typename MatrixType::type_t Scalar; 379 typedef typename NumTraits<Scalar>::ComplexScalar ComplexScalar; 380 381 typedef VectorEigenDense<ComplexScalar> EigenvalueType; 382 typedef VectorEigenDense<Scalar> ColumnVectorType; 383 384 /** \brief Default constructor. 385 * 386 * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. 387 * 388 * The default constructor is useful in cases in which the user intends to 389 * perform decompositions via compute(). The \p size parameter is only 390 * used as a hint. It is not an error to give a wrong \p size, but it may 391 * impair performance. 392 * 393 * \sa compute() for an example. 394 */ 395 RealSchur(Index size) 396 : matT_(size, size), 397 matU_(size, size), 398 hess_(size), 399 maxIterations_(40), 400 isInitialized_(false), 401 matUisUptodate_(false), 402 info_(_noConvergence) 403 { } 404 ... 405\end{lstlisting} 406\vspace{.1cm} 407We can use the constructor which computes the real Schur decomposition at construction time. 408\vspace{.1cm} 409\begin{lstlisting}[]{} 410RealSchur(const MatrixType& matrix, bool computeU = true); 411\end{lstlisting} 412\vspace{.1cm} 413Once the decomposition is computed, we can use the {\cmd matrixU()} and {\cmd matrixT()} functions to retrieve the matrices $U$ and $T$ in the decomposition. 414\vspace{.1cm} 415\begin{lstlisting}[]{} 416MatrixType matrixU() const; 417MatrixType matrixT() const; 418\end{lstlisting} 419\vspace{.1cm} 420Alternatively calling the function {\cmd compute()} to compute the real Schur decomposition of a given matrix, after that we are able to use these two functions above. \\ 421The implementation of Schur decomposition of real square matrix is adapted from \url{http://math.nist.gov/javanumerics/jama/}, whose code is based on \eispack. The Schur decomposition is computed by first reducing the matrix to $Hessenberg$ form using the class {\class HessenbergDecomposition}. The $Hessenberg$ matrix is then reduced to triangular form by performing $Francis QR$ iterations with implicit double shift.\\ 422There are several auxiliary functions in order to implement this algorithm : 423\vspace{.1cm} 424\begin{lstlisting}[deletekeywords={[3] norm}] 425Scalar computeNormOfT(); 426Index findSmallSubdiagEntry(Index iu, Scalar norm); 427void splitOffTwoRows(Index iu, bool computeU, Scalar exshift); 428void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); 429void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); 430void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector); 431\end{lstlisting} 432\vspace{.1cm} 433For instance, the following $5\times 5$ non symmetric matrix 434\[ 435A=\left[ 436\begin{array}{cccccc} 437 1.8779& -0.5583& -0.2099& 0.2696& 0.1097\\ 4380.9407& -0.3114& -1.6989& 0.4943& 1.1287 \\ 4390.7873& -0.57& 0.6076& -1.4831& -0.29 \\ 440-0.8759& -1.0257& -0.1178& -1.0203& 1.2616 \\ 4410.3199& -0.9087& 0.6992& -0.447& 0.4754 442\end{array} 443\right] 444\] 445has the quasi-triangular matrix in the Schur decomposition 446\footnotesize 447\[ 448H=\left[ 449\begin{array}{cccccc} 4501.6552531332& -0.30767059895 & -0.45435806468 & 1.5602779282 & 0.68358491675\\ 451 0 & 1.1891303949 & 1.9306398993 & 0.96785347442 & 0.57970343469 \\ 452 0 & -1.3969331077& 0.22387492903 & -0.78380738142& 1.27402357 \\ 453 0 & 0 & 0 & 0.092450608769& -0.5563858236 \\ 454 0 & 0 & 0 & 0 & -1.5315090659 455\end{array} 456\right] 457\] 458\normalsize 459and the orthogonal matrix in the Schur decomposition. 460\footnotesize 461\[ 462Q=\left[ 463\begin{array}{cccccc} 4640.72665944996& 0.28908942548& -0.090578168382& 0.61658006012& 0.0042394246225 \\ 465 0.13365823926& 0.75252664471 & -0.11720311204 & -0.52512343394& -0.35545454525 \\ 466 0.50462706182 & -0.59105620177& -0.20754951129& -0.34476029525& -0.4838053679 \\ 467 -0.1068784212 & 0.010284495663 & 0.747282602 & 0.23512571367& -0.61217305049 \\ 468 0.43362431945& -0.026033604674& 0.61363746549 & -0.41222645012& 0.51459752796 469\end{array} 470\right] 471\] 472\normalsize 473\displayInfos{library=eigenSolvers, test={test\_EigenSolverDecomposition.cpp}, header dep={Jacobi.hpp, VectorEigenDense.hpp, HessenbergDecomposition.hpp}} 474 475\subsection{{\libtitle eigenSolver} main part} 476 477In this section, there are solvers for three kinds of eigen problem : 478\begin{itemize} 479\item Hermitian eigenproblem 480\item Generalized Hermitian eigenproblem 481\item Non-Hermitian (real) eigenproblem 482\end{itemize} 483The following presents eigen solver classes corresponding to each kind of eigen problem. 484 485\subsubsection{The {\classtitle SelfAdjointEigenSolver} class} 486 487Class {\class SelfAdjointEigenSolver} computes eigenvalues and eigenvectors of self-adjoint matrices. 488A matrix $A$ is selfadjoint if it equals its adjoint. For real matrices, this means that the matrix is symmetric: it equals its transpose. \\ 489This class computes the eigenvalues and eigenvectors of a self-adjoint matrix. These are the scalars $ \lambda $ and vectors $ v $ such that $ Av = \lambda v $. The eigenvalues of a self-adjoint matrix are always real. If $ D $ is a diagonal matrix with the eigenvalues on the diagonal, and $ V $ is a matrix with the eigenvectors as its columns, then $ A = V D V^{-1} $ (for self-adjoint matrices, the matrix $ V $ is always invertible). This is called the eigendecomposition. The implemented algorithm exploits the fact that the matrix is self-adjoint, making it faster and more accurate than the general purpose eigenvalue algorithms implemented in RealEigenSolver.\\ 490\vspace{.1cm} 491\begin{lstlisting} 492template<typename _MatrixType> 493class SelfAdjointEigenSolver 494{ 495 public: 496 497 typedef _MatrixType MatrixType; 498 499 /** \brief Scalar type for matrices of type \p _MatrixType. */ 500 typedef typename MatrixType::type_t Scalar; 501 502 /** \brief Real scalar type for \p _MatrixType. 503 * 504 * This is just \c Scalar if #Scalar is real (e.g., \c float or 505 * \c double), and the type of the real part of \c Scalar if #Scalar is 506 * complex. 507 */ 508 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 509 510 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 511 * 512 * This is a column vector with entries of type #RealScalar. 513 * The length of the vector is the size of \p _MatrixType. 514 */ 515 typedef VectorEigenDense<RealScalar> RealVectorType; 516 typedef Tridiagonalization<MatrixType> TridiagonalizationType; 517 518 /** \brief Default constructor for fixed-size matrices. 519 * 520 * The default constructor is useful in cases in which the user intends to 521 * perform decompositions via compute(). 522 * 523 */ 524 SelfAdjointEigenSolver() 525 : eivec_(), 526 eivalues_(), 527 subdiag_(), 528 maxIterations_(30), 529 isInitialized_(false) 530 { } 531 ... 532\end{lstlisting} 533\vspace{.1cm} 534We can use the constructor which computes the eigenvalues and eigenvectors at construction time. 535\vspace{.1cm} 536\begin{lstlisting}[]{} 537SelfAdjointEigenSolver(const MatrixType& matrix, Int options = _computeEigenVector); 538\end{lstlisting} 539\vspace{.1cm} 540Once the eigenvalue and eigenvectors are computed, they can be retrieved with 541\vspace{.1cm} 542\begin{lstlisting}[]{} 543const RealVectorType& eigenvalues() const; 544const MatrixType& eigenvectors() const; 545\end{lstlisting} 546\vspace{.1cm} 547Alternatively, we can call the function {\cmd compute()} to compute the eigenvalues and eigenvectors of a given matrix. 548\vspace{.1cm} 549\begin{lstlisting}[]{} 550SelfAdjointEigenSolver& compute(const MatrixType& matrix, Int options = _computeEigenVector); 551\end{lstlisting} 552\vspace{.1cm} 553This implementation uses a symmetric QR algorithm. The matrix is first reduced to tridiagonal form using the {\class Tridiagonalization} class. The tridiagonal matrix is then brought to diagonal form with implicit symmetric QR steps with Wilkinson shift.\\ 554For instance, the following $5\times 5$ real symmetric matrix 555\[ 556A=\left[ 557\begin{array}{cccccc} 5581.36 & -0.816 & 0.521 & 1.43 & -0.144\\ 559-0.816 & -0.659 & 0.794 & -0.173 & -0.406 \\ 5600.521 & 0.794 & -0.541 & 0.461& 0.179 \\ 5611.43& -0.173 & 0.461& -1.43 & 0.822 \\ 562-0.144 & -0.406 & 0.179 & 0.822& -1.37 563\end{array} 564\right] 565\] 566has eigenvalues and corresponding (column) eigenvector\\ 567\mbox{}\hspace{10 mm} eigenvalues = ( -2.64925 -1.76682 -0.742793 0.22723 2.29163 ) \\ 568\footnotesize 569\[ 570\mbox{\normalsize eigenvectors =} \left[ 571\begin{array}{cccccc} 572-0.32648814917& -0.098193931162& 0.34709885718& -0.010999922777& 0.87359305479 \\ 573-0.20806598152& -0.64229394812& 0.22733262857 & 0.66233712946 & -0.23194058392 \\ 5740.050094522003 & 0.62960023585& -0.16376394055 & 0.73987891396& 0.16387387324 \\ 5750.72069515152 & -0.39678270513 & -0.40156098976& 0.11452959441& 0.38573789934 \\ 576-0.57289010572& -0.15486595502& -0.79985773898& -0.025508690809 & 0.085967242812 577\end{array} 578\right] 579\] 580\normalsize 581 582\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={VectorEigenDense.hpp, Tridiagonalization.hpp}} 583 584\subsubsection{The {\classtitle GeneralizedSelfAdjointEigenSolver} class} 585 586Class {\class GeneralizedSelfAdjointEigenSolver} computes eigenvalues and eigenvectors of the generalized self-adjoint eigen problem. \\ This class solves the generalized eigenvalue problem $ Av = \lambda Bv $. In this case, the matrix $ A $ should be self-adjoint and the matrix $ B $ should be positive definite. 587\vspace{.1cm} 588\begin{lstlisting}[deletekeywords={[3] size}] 589template<typename _MatrixType> 590class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> 591{ 592 typedef SelfAdjointEigenSolver<_MatrixType> Base; 593public: 594 typedef _MatrixType MatrixType; 595 596 /** \brief Default constructor for fixed-size matrices. 597 * 598 * The default constructor is useful in cases in which the user intends to 599 * perform decompositions via compute(). 600 */ 601 GeneralizedSelfAdjointEigenSolver() : Base() {} 602 603 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 604 * 605 * \param [in] size Positive integer, size of the matrix whose 606 * eigenvalues and eigenvectors will be computed. 607 * 608 * \sa compute() for an example 609 */ 610 GeneralizedSelfAdjointEigenSolver(Index size) 611 : Base(size) 612 {} 613 ... 614\end{lstlisting} 615\vspace{.1cm} 616We can use the constructor which computes the eigenvalues and eigenvectors at construction time. 617\vspace{.1cm} 618\begin{lstlisting}[]{} 619 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 620 int options = _computeEigenVector|_Ax_lBx) 621 : Base(matA.numOfCols()); 622\end{lstlisting} 623\vspace{.1cm} 624Alternatively, we can call the function {\itshape compute()} to compute the eigenvalues and eigenvectors of a given matrix. 625\vspace{.1cm} 626\begin{lstlisting}[]{} 627 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 628 int options = _computeEigenVector|_Ax_lBx); 629\end{lstlisting} 630\vspace{.1cm} 631Eigenvalues and eigenvectors can be retrieved with same functions of {\class SelfAdjointEigenSolver} 632\vspace{.1cm} 633\begin{lstlisting}[]{} 634const RealVectorType& eigenvalues() const; 635const MatrixType& eigenvectors() const; 636\end{lstlisting} 637\vspace{.1cm} 638For instance, the following $5\times 5$ real symmetric matrix $A$ and positive definite matrix $B$ 639\[ 640A=\left[ 641\begin{array}{cccccc} 6421.36 & -0.816 & 0.521 & 1.43 & -0.144\\ 643-0.816 & -0.659 & 0.794 & -0.173 & -0.406 \\ 6440.521 & 0.794 & -0.541 & 0.461& 0.179 \\ 6451.43& -0.173 & 0.461& -1.43 & 0.822 \\ 646-0.144 & -0.406 & 0.179 & 0.822& -1.37 647\end{array} 648\right] 649\] 650\[ 651B=\left[ 652\begin{array}{cccccc} 6531.1981& 1.4405& 1.5367& 1.0971& 1.1276\\ 6541.4405& 2.1207& 1.9115& 1.3836& 1.642\\ 6551.5367& 1.9115& 2.4828& 1.7702& 1.7543\\ 6561.0971& 1.3836& 1.7702& 1.4359& 1.2428\\ 6571.1276& 1.642& 1.7543& 1.2428& 1.424 658\end{array} 659\right] 660\] 661has eigenvalues and corresponding (column) eigenvector\\ 662\mbox{}\hspace{10 mm} eigenvalues = ( -40.0667 -12.6587 -3.08602 0.0522151 40.3448 ) \\ 663\footnotesize 664\[ 665\mbox{\normalsize eigenvectors =} \left[ 666\begin{array}{cccccc} 667-2.5241477328& -0.01326907506 & 0.83333293087 & -0.040348533316 & 4.3939678718\\ 6683.7607674197 & 0.3488243247 & 1.1073133113 & -0.3174355818 & -4.3303408407 \\ 6693.0396247156 & 2.615639864 & -0.34676637691 & -0.36756457169 & -3.9258019948 \\ 6700.45843793898 & -2.0051996759 & -0.65863554474 & -0.039469282032 & 1.4762595546 \\ 671-6.6376554636 & -1.9381010232 & -0.97022016252 & 0.070809133211 & 4.9718019236 672\end{array} 673\right] 674\] 675\normalsize 676\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={Tridiagonalization.hpp}} 677 678\subsubsection{The {\classtitle RealEigenSolver} class} 679 680Class {\class RealEigenSolver} computes eigenvalues and eigenvectors of general real matrices. \\ 681The eigenvalues and eigenvectors of a matrix $ A $ are scalars $ \lambda $ and vectors $ v $ such that $ Av = \lambda v $. If $ D $ is a diagonal matrix with the eigenvalues on the diagonal, and $ V $ is a matrix with the eigenvectors as its columns, then $ A V = V D $. The matrix $ V $ is almost always invertible, in which case we have $ A = V D V^{-1} $. \\ 682The eigenvalues and eigenvectors of the matrix may be complex, even when the matrix is real. However, we can choose real matrices $ V $ and $ D $ satisfying $ A V = V D $, just like the eigendecomposition, if the matrix $ D $ is not required to be diagonal, but if it is allowed to have blocks of the form\[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \] (where $ u $ and $ v $ are real numbers) on the diagonal. These blocks correspond to complex eigenvalue pairs $ u \pm iv $. We call this variant of the eigendecomposition the pseudo-eigendecomposition. 683\vspace{.1cm} 684\begin{lstlisting}[deletekeywords={[3] size}] 685template<typename _MatrixType> 686class RealEigenSolver 687{ 688 public: 689 690 /** \brief Synonym for the template parameter \p _MatrixType. */ 691 typedef _MatrixType MatrixType; 692 693 /** \brief Scalar type for matrices of type #MatrixType. */ 694 typedef typename MatrixType::type_t Scalar; 695 typedef typename NumTraits<Scalar>::RealScalar RealScalar; 696 typedef VectorEigenDense<Scalar> ColumnVectorType; 697 typedef VectorEigenDense<Scalar> VectorType; 698 699 /** \brief Complex scalar type for #MatrixType. 700 * 701 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 702 * \c float or \c double) and just \c Scalar if #Scalar is 703 * complex. 704 */ 705 typedef typename NumTraits<Scalar>::ComplexScalar ComplexScalar; 706 typedef VectorEigenDense<ComplexScalar> ComplexVectorType; 707 708 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 709 * 710 * This is a column vector with entries of type #ComplexScalar. 711 * The length of the vector is the size of #MatrixType. 712 */ 713 typedef VectorEigenDense<ComplexScalar> EigenvalueType; 714 715 /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 716 * 717 * This is a square matrix with entries of type #ComplexScalar. 718 * The size is the same as the size of #MatrixType. 719 */ 720 typedef MatrixEigenDense<ComplexScalar> EigenvectorsType; 721 722 /** \brief Default constructor. 723 * 724 * The default constructor is useful in cases in which the user intends to 725 * perform decompositions via RealEigenSolver::compute(const MatrixType&, bool). 726 * 727 * \sa compute() for an example. 728 */ 729 RealEigenSolver() : eivec_(), eivalues_(), isInitialized_(false), realSchur_(), matT_() {} 730 731 /** 732 * \brief Default constructor with memory preallocation 733 * 734 * Like the default constructor but with preallocation of the internal data 735 * according to the specified problem \a size. 736 * \sa RealEigenSolver() 737 */ 738 RealEigenSolver(Index size) 739 : eivec_(size, size), 740 eivalues_(size), 741 isInitialized_(false), 742 eigenvectorsOk_(false), 743 realSchur_(size), 744 matT_(size, size) 745 {} 746 ... 747\end{lstlisting} 748\vspace{.1cm} 749Same as other two eigensolvers, we can use the constructor which computes the eigenvalues and eigenvectors at construction time. 750\vspace{.1cm} 751\begin{lstlisting}[]{} 752RealEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true); 753\end{lstlisting} 754\vspace{.1cm} 755Once the eigenvalue and eigenvectors are computed, they can be retrieved with 756\vspace{.1cm} 757\begin{lstlisting}[]{} 758const EigenvalueType& eigenvalues() const; 759EigenvectorsType eigenvectors() const; 760\end{lstlisting} 761\vspace{.1cm} 762Alternatively, we can call the function {\cmd compute()} to compute the eigenvalues and eigenvectors of a given matrix. 763\vspace{.1cm} 764\begin{lstlisting}[]{} 765RealEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 766\end{lstlisting} 767\vspace{.1cm} 768This implementation first reduces the matrix to real Schur form using the {\class RealSchur} class. The Schur decomposition is then used to compute the eigenvalues and eigenvectors.\\ 769For instance, the following $5\times 5$ non symmetric matrix 770\[ 771A=\left[ 772\begin{array}{cccccc} 7731.8779& -0.5583& -0.2099& 0.2696& 0.1097\\ 7740.9407& -0.3114& -1.6989& 0.4943& 1.1287 \\ 7750.7873& -0.57& 0.6076& -1.4831& -0.29 \\ 776-0.8759& -1.0257& -0.1178& -1.0203& 1.2616 \\ 7770.3199& -0.9087& 0.6992& -0.447& 0.4754 778\end{array} 779\right] 780\] 781has eigenvalues and corresponding (column) eigenvector\\ 782\mbox{}\hspace{10 mm} eigenvalues = ( (1.65525,0) (0.706503,1.56973) (0.706503,-1.56973) (0.0924506,0) (-1.53151,0) ) \\ 783\[ 784\mbox{eigenvectors =} \left[ 785\begin{array}{ccccc} 786(0.727,0)& (0.177,-0.001)& (0.177,0.001)& (-0.230,0)& (-0.041,0) \\ 787(0.133,0)& (0.534,-0.218)& (0.534,0.218)& (-0.664,0)& (-0.392,0) \\ 788(0.504,0)& (-0.444,0.0814) & (-0.444,-0.0814)& (-0.336,0)& (-0.587,0)\\ 789(-0.106,0)& (0.011,0.457)& (0.0118,-0.457)& (0.135,0)& (-0.695,0) \\ 790(0.433,0)& (-0.036,0.469)& (-0.036,-0.469)& (-0.611,0)& (-0.121,0) 791\end{array} 792\right] 793\] 794 795\displayInfos{library=eigenSolvers, test={test\_EigenSolverInternCore.cpp}, header dep={VectorEigenDense.hpp, MatrixEigenDense.hpp, RealSchur.hpp}} 796