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