Michael Ernst Peter / Eigen

src/SVD/UpperBidiagonalization.h@1:3b8049da21b8, 2019-09-24 (annotated)

Committer:
jsoh91
Date:
Tue Sep 24 00:18:23 2019 +0000
Revision:
1:3b8049da21b8
Parent:
0:13a5d365ba16 
ignore and revise some of error parts

Who changed what in which revision?

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