UpperBidiagonalization.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_BIDIAGONALIZATION_H
26 #define EIGEN_BIDIAGONALIZATION_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 // UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
32 // At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.
33 
34 template<typename _MatrixType> class UpperBidiagonalization
35 {
36  public:
37 
38  typedef _MatrixType MatrixType;
39  enum {
40  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
41  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
42  ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
43  };
44  typedef typename MatrixType::Scalar Scalar;
45  typedef typename MatrixType::RealScalar RealScalar;
46  typedef typename MatrixType::Index Index;
47  typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
48  typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
49  typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType;
50  typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
51  typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
52  typedef HouseholderSequence<
53  const MatrixType,
54  CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, const Diagonal<const MatrixType,0> >
55  > HouseholderUSequenceType;
56  typedef HouseholderSequence<
57  const MatrixType,
58  Diagonal<const MatrixType,1>,
60  > HouseholderVSequenceType;
61 
68  UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
69 
70  UpperBidiagonalization(const MatrixType& matrix)
71  : m_householder(matrix.rows(), matrix.cols()),
72  m_bidiagonal(matrix.cols(), matrix.cols()),
73  m_isInitialized(false)
74  {
75  compute(matrix);
76  }
77 
78  UpperBidiagonalization& compute(const MatrixType& matrix);
79 
80  const MatrixType& householder() const { return m_householder; }
81  const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
82 
83  const HouseholderUSequenceType householderU() const
84  {
85  eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
86  return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
87  }
88 
89  const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
90  {
91  eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
92  return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
93  .setLength(m_householder.cols()-1)
94  .setShift(1);
95  }
96 
97  protected:
98  MatrixType m_householder;
99  BidiagonalType m_bidiagonal;
100  bool m_isInitialized;
101 };
102 
103 template<typename _MatrixType>
104 UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
105 {
106  Index rows = matrix.rows();
107  Index cols = matrix.cols();
108 
109  eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols.");
110 
111  m_householder = matrix;
112 
113  ColVectorType temp(rows);
114 
115  for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
116  {
117  Index remainingRows = rows - k;
118  Index remainingCols = cols - k - 1;
119 
120  // construct left householder transform in-place in m_householder
121  m_householder.col(k).tail(remainingRows)
122  .makeHouseholderInPlace(m_householder.coeffRef(k,k),
123  m_bidiagonal.template diagonal<0>().coeffRef(k));
124  // apply householder transform to remaining part of m_householder on the left
125  m_householder.bottomRightCorner(remainingRows, remainingCols)
126  .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1),
127  m_householder.coeff(k,k),
128  temp.data());
129 
130  if(k == cols-1) break;
131 
132  // construct right householder transform in-place in m_householder
133  m_householder.row(k).tail(remainingCols)
134  .makeHouseholderInPlace(m_householder.coeffRef(k,k+1),
135  m_bidiagonal.template diagonal<1>().coeffRef(k));
136  // apply householder transform to remaining part of m_householder on the left
137  m_householder.bottomRightCorner(remainingRows-1, remainingCols)
138  .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(),
139  m_householder.coeff(k,k+1),
140  temp.data());
141  }
142  m_isInitialized = true;
143  return *this;
144 }
145 
146 #if 0
147 
151 template<typename Derived>
152 const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
153 MatrixBase<Derived>::bidiagonalization() const
154 {
155  return UpperBidiagonalization<PlainObject>(eval());
156 }
157 #endif
158 
159 } // end namespace internal
160 
161 } // end namespace Eigen
162 
163 #endif // EIGEN_BIDIAGONALIZATION_H