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
13namespace Eigen {
14
15namespace 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
19template<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    */
53    UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}
54
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
65    const MatrixType& householder() const { return m_householder; }
66    const BidiagonalType& bidiagonal() const { return m_bidiagonal; }
67
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
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
88template<typename _MatrixType>
89UpperBidiagonalization<_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  */
136template<typename Derived>
137const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
138MatrixBase<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