1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_TRIDIAGONALIZATION_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_TRIDIAGONALIZATION_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct TridiagonalizationMatrixTReturnType; 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<TridiagonalizationMatrixTReturnType<MatrixType> > 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public traits<typename MatrixType::PlainObject> 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix? 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { Flags = 0 }; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename CoeffVectorType> 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs); 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class Tridiagonalization 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Tridiagonal decomposition of a selfadjoint matrix 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are computing the 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * tridiagonal decomposition; this is expected to be an instantiation of the 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Matrix class template. 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that: 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix. 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A tridiagonal matrix is a matrix which has nonzero elements only on the 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * main diagonal and the first diagonal below and above it. The Hessenberg 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition of a selfadjoint matrix is in fact a tridiagonal 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition. This class is used in SelfAdjointEigenSolver to compute the 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues and eigenvectors of a selfadjoint matrix. 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Call the function compute() to compute the tridiagonal decomposition of a 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&) 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * constructor which computes the tridiagonal Schur decomposition at 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * construction time. Once the decomposition is computed, you can use the 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition. 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The documentation of Tridiagonalization(const MatrixType&) contains an 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * example of the typical use of this class. 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class Tridiagonalization 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Synonym for the template parameter \p _MatrixType. */ 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Size = MatrixType::RowsAtCompileTime, 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1), 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxSize = MatrixType::MaxRowsAtCompileTime, 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1) 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type, 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Diagonal<const MatrixType> 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath >::type DiagonalReturnType; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type, 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Diagonal<const MatrixType, -1> 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath >::type SubDiagonalReturnType; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Return type of matrixQ() */ 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType; 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor. 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param [in] size Positive integer, size of the matrix whose tridiagonal 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition will be computed. 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via compute(). The \p size parameter is only 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * used as a hint. It is not an error to give a wrong \p size, but it may 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * impair performance. 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example. 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size) 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(size,size), 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs(size > 1 ? size-1 : 1), 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes tridiagonal decomposition of given matrix. 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is to be computed. 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute() to compute the tridiagonal decomposition. 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit Tridiagonalization(const EigenBase<InputType>& matrix) 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_matrix(matrix.derived()), 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1), 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::tridiagonalization_inplace(m_matrix, m_hCoeffs); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes tridiagonal decomposition of given matrix. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is to be computed. 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The tridiagonal decomposition is computed by bringing the columns of 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the matrix successively in the required form using Householder 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the size of the given matrix. 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method reuses of the allocated data in the Tridiagonalization 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * object, if the size of the matrix does not change. 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Tridiagonalization_compute.cpp 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_compute.out 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Tridiagonalization& compute(const EigenBase<InputType>& matrix) 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matrix = matrix.derived(); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs.resize(matrix.rows()-1, 1); 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::tridiagonalization_inplace(m_matrix, m_hCoeffs); 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the Householder coefficients. 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a const reference to the vector of Householder coefficients 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Householder coefficients allow the reconstruction of the matrix 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q \f$ in the tridiagonal decomposition from the packed data. 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Tridiagonalization_householderCoefficients.cpp 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_householderCoefficients.out 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa packedMatrix(), \ref Householder_Module "Householder module" 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline CoeffVectorType householderCoefficients() const 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_hCoeffs; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the internal representation of the decomposition 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a const reference to a matrix with the internal representation 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the decomposition. 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The returned matrix contains the following information: 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - the strict upper triangular part is equal to the input matrix A. 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - the diagonal and lower sub-diagonal represent the real tridiagonal 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * symmetric matrix T. 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - the rest of the lower part contains the Householder vectors that, 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * combined with Householder coefficients returned by 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * householderCoefficients(), allows to reconstruct the matrix Q as 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Here, the matrices \f$ H_i \f$ are the Householder transformations 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H_i = (I - h_i v_i v_i^T) \f$ 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i \f$ is the Householder vector defined by 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$ 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with M the matrix returned by this function. 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See LAPACK for further details on this packed storage. 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Tridiagonalization_packedMatrix.cpp 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_packedMatrix.out 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa householderCoefficients() 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixType& packedMatrix() const 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrix; 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the unitary matrix Q in the decomposition 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns object representing the matrix Q 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a light-weight object of template class 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HouseholderSequence. You can either apply it directly to a matrix or 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * you can convert it to a matrix of type #MatrixType. 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa Tridiagonalization(const MatrixType&) for an example, 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrixT(), class HouseholderSequence 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HouseholderSequenceType matrixQ() const 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setLength(m_matrix.rows() - 1) 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setShift(1); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns an expression of the tridiagonal matrix T in the decomposition 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns expression object representing the matrix T 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Currently, this function can be used to extract the matrix T from internal 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * data and copy it to a dense matrix object. In most cases, it may be 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * sufficient to directly use the packed matrix or the vector expressions 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * returned by diagonal() and subDiagonal() instead of creating a new 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * dense copy matrix with this function. 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa Tridiagonalization(const MatrixType&) for an example, 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrixQ(), packedMatrix(), diagonal(), subDiagonal() 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixTReturnType matrixT() const 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return MatrixTReturnType(m_matrix.real()); 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition. 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns expression representing the diagonal of T 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Tridiagonalization_diagonal.cpp 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_diagonal.out 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa matrixT(), subDiagonal() 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagonalReturnType diagonal() const; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition. 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns expression representing the subdiagonal of T 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor Tridiagonalization(const MatrixType&) or 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the member function compute(const MatrixType&) has been called before 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the tridiagonal decomposition of a matrix. 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa diagonal() for an example, matrixT() 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SubDiagonalReturnType subDiagonal() const; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matrix; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CoeffVectorType m_hCoeffs; 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename Tridiagonalization<MatrixType>::DiagonalReturnType 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathTridiagonalization<MatrixType>::diagonal() const 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 3102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_matrix.diagonal().real(); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename Tridiagonalization<MatrixType>::SubDiagonalReturnType 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathTridiagonalization<MatrixType>::subDiagonal() const 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_matrix.template diagonal<-1>().real(); 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place. 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced. 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On output, the strict upper part is left unchanged, and the lower triangular part 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * represents the T and Q matrices in packed format has detailed below. 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] hCoeffs returned Householder coefficients (see below) 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On output, the tridiagonal selfadjoint matrix T is stored in the diagonal 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and lower sub-diagonal of the matrix \a matA. 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The unitary matrix Q is represented in a compact way as a product of 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Householder reflectors \f$ H_i \f$ such that: 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Householder reflectors are defined as 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H_i = (I - h_i v_i v_i^T) \f$ 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where \f$ h_i = hCoeffs[i]\f$ is the \f$ i \f$th Householder coefficient and 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i \f$ is the Householder vector defined by 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$. 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa Tridiagonalization::packedMatrix() 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename CoeffVectorType> 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using numext::conj; 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = matA.rows(); 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(n==matA.cols()); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(n==hCoeffs.size()+1 || n==1); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i<n-1; ++i) 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index remainingSize = n-i-1; 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar beta; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar h; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply similarity transformation to remaining columns, 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1) 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(i).coeffRef(i+1) = 1; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>() 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (conj(h) * matA.col(i).tail(remainingSize))); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang hCoeffs.tail(n-i-1) += (conj(h)*RealScalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>() 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1)); 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(i).coeffRef(i+1) = beta; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hCoeffs.coeffRef(i) = h; 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// forward declaration, implementation at the end of this file 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int Size=MatrixType::ColsAtCompileTime, 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex> 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct tridiagonalization_inplace_selector; 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \brief Performs a full tridiagonalization in place 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in,out] mat On input, the selfadjoint matrix whose tridiagonal 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition is to be computed. Only the lower triangular part referenced. 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The rest is left unchanged. On output, the orthogonal matrix Q 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * in the decomposition if \p extractQ is true. 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] diag The diagonal of the tridiagonal matrix T in the 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition. 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] subdiag The subdiagonal of the tridiagonal matrix T in 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the decomposition. 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] extractQ If true, the orthogonal matrix Q in the 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition is computed and stored in \p mat. 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes the tridiagonal decomposition of the selfadjoint matrix \p mat in place 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * such that \f$ mat = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * symmetric tridiagonal matrix. 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The tridiagonal matrix T is passed to the output parameters \p diag and \p subdiag. If 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat. Otherwise the lower 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part of the matrix \p mat is destroyed. 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors \p diag and \p subdiag are not resized. The function 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * assumes that they are already of the correct size. The length of the 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * vector \p diag should equal the number of rows in \p mat, and the 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * length of the vector \p subdiag should be one left. 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This implementation contains an optimized path for 3-by-3 matrices 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * which is especially useful for plane fitting. 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note Currently, it requires two temporary vectors to hold the intermediate 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Householder coefficients, and to reconstruct the matrix Q from the Householder 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * reflectors. 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example (this uses the same matrix as the example in 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Tridiagonalization::Tridiagonalization(const MatrixType&)): 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \include Tridiagonalization_decomposeInPlace.cpp 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Tridiagonalization_decomposeInPlace.out 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class Tridiagonalization 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename DiagonalType, typename SubDiagonalType> 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * General full tridiagonalization 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int Size, bool IsComplex> 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct tridiagonalization_inplace_selector 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType; 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType; 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DiagonalType, typename SubDiagonalType> 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CoeffVectorType hCoeffs(mat.cols()-1); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tridiagonalization_inplace(mat,hCoeffs); 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag = mat.diagonal().real(); 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath subdiag = mat.template diagonal<-1>().real(); 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(extractQ) 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat = HouseholderSequenceType(mat, hCoeffs.conjugate()) 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setLength(mat.rows() - 1) 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setShift(1); 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Specialization for 3x3 real matrices. 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Especially useful for plane fitting. 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct tridiagonalization_inplace_selector<MatrixType,3,false> 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DiagonalType, typename SubDiagonalType> 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 4692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[0] = mat(0,0); 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar v1norm2 = numext::abs2(mat(2,0)); 4722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(v1norm2 <= tol) 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[1] = mat(1,1); 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[2] = mat(2,2); 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath subdiag[0] = mat(1,0); 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath subdiag[1] = mat(2,1); 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (extractQ) 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.setIdentity(); 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2); 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar invBeta = RealScalar(1)/beta; 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar m01 = mat(1,0) * invBeta; 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar m02 = mat(2,0) * invBeta; 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1)); 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[1] = mat(1,1) + m02*q; 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath diag[2] = mat(2,2) - m02*q; 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath subdiag[0] = beta; 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath subdiag[1] = mat(2,1) - m01 * q; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (extractQ) 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat << 1, 0, 0, 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, m01, m02, 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, m02, -m01; 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Trivial specialization for 1x1 matrices 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, bool IsComplex> 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct tridiagonalization_inplace_selector<MatrixType,1,IsComplex> 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename DiagonalType, typename SubDiagonalType> 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ) 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez diag(0,0) = numext::real(mat(0,0)); 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(extractQ) 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat(0,0) = Scalar(1); 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \eigenvalues_module \ingroup Eigenvalues_Module 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Expression type for return value of Tridiagonalization::matrixT() 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of underlying dense matrix 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct TridiagonalizationMatrixTReturnType 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath: public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> > 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] mat The underlying dense matrix 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { } 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void evalTo(ResultType& result) const 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.setZero(); 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate(); 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.diagonal() = m_matrix.diagonal(); 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.template diagonal<-1>() = m_matrix.template diagonal<-1>(); 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_matrix.rows(); } 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_matrix.cols(); } 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Nested m_matrix; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_TRIDIAGONALIZATION_H 557