1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_BIDIAGONALIZATION_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_BIDIAGONALIZATION_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API. 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class. 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class UpperBidiagonalization 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef HouseholderSequence< 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType, 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> > 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath > HouseholderUSequenceType; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef HouseholderSequence< 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type, 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Diagonal<const MatrixType,1>, 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath OnTheRight 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath > HouseholderVSequenceType; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Default Constructor. 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via Bidiagonalization::compute(const MatrixType&). 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {} 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpperBidiagonalization(const MatrixType& matrix) 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_householder(matrix.rows(), matrix.cols()), 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_bidiagonal(matrix.cols(), matrix.cols()), 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpperBidiagonalization& compute(const MatrixType& matrix); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& householder() const { return m_householder; } 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const BidiagonalType& bidiagonal() const { return m_bidiagonal; } 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const HouseholderUSequenceType householderU() const 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate()); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized."); 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>()) 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setLength(m_householder.cols()-1) 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setShift(1); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_householder; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BidiagonalType m_bidiagonal; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathUpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = matrix.rows(); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = matrix.cols(); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols."); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder = matrix; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColVectorType temp(rows); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k) 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index remainingRows = rows - k; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index remainingCols = cols - k - 1; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // construct left householder transform in-place in m_householder 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.col(k).tail(remainingRows) 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .makeHouseholderInPlace(m_householder.coeffRef(k,k), 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_bidiagonal.template diagonal<0>().coeffRef(k)); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // apply householder transform to remaining part of m_householder on the left 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.bottomRightCorner(remainingRows, remainingCols) 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1), 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.coeff(k,k), 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp.data()); 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k == cols-1) break; 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // construct right householder transform in-place in m_householder 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.row(k).tail(remainingCols) 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .makeHouseholderInPlace(m_householder.coeffRef(k,k+1), 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_bidiagonal.template diagonal<1>().coeffRef(k)); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // apply householder transform to remaining part of m_householder on the left 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.bottomRightCorner(remainingRows-1, remainingCols) 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(), 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_householder.coeff(k,k+1), 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp.data()); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#if 0 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \return the Householder QR decomposition of \c *this. 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class Bidiagonalization 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject> 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::bidiagonalization() const 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return UpperBidiagonalization<PlainObject>(eval()); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_BIDIAGONALIZATION_H 149