Another set of UI messages fixes and tweaks! No functional changes.
[blender.git] / extern / Eigen3 / Eigen / src / Householder / HouseholderSequence.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25
26 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
27 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
28
29 /** \ingroup Householder_Module
30   * \householder_module
31   * \class HouseholderSequence
32   * \brief Sequence of Householder reflections acting on subspaces with decreasing size
33   * \tparam VectorsType type of matrix containing the Householder vectors
34   * \tparam CoeffsType  type of vector containing the Householder coefficients
35   * \tparam Side        either OnTheLeft (the default) or OnTheRight
36   *
37   * This class represents a product sequence of Householder reflections where the first Householder reflection
38   * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by
39   * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace
40   * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but
41   * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections
42   * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods
43   * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(),
44   * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence.
45   *
46   * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the
47   * form \f$ H = \prod_{i=0}^{n-1} H_i \f$ where the i-th Householder reflection is \f$ H_i = I - h_i v_i
48   * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$
49   * v_i \f$ is a vector of the form
50   * \f[ 
51   * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
52   * \f]
53   * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector.
54   *
55   * Typical usages are listed below, where H is a HouseholderSequence:
56   * \code
57   * A.applyOnTheRight(H);             // A = A * H
58   * A.applyOnTheLeft(H);              // A = H * A
59   * A.applyOnTheRight(H.adjoint());   // A = A * H^*
60   * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
61   * MatrixXd Q = H;                   // conversion to a dense matrix
62   * \endcode
63   * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators.
64   *
65   * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example.
66   *
67   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
68   */
69
70 namespace internal {
71
72 template<typename VectorsType, typename CoeffsType, int Side>
73 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
74 {
75   typedef typename VectorsType::Scalar Scalar;
76   typedef typename VectorsType::Index Index;
77   typedef typename VectorsType::StorageKind StorageKind;
78   enum {
79     RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
80                                         : traits<VectorsType>::ColsAtCompileTime,
81     ColsAtCompileTime = RowsAtCompileTime,
82     MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
83                                            : traits<VectorsType>::MaxColsAtCompileTime,
84     MaxColsAtCompileTime = MaxRowsAtCompileTime,
85     Flags = 0
86   };
87 };
88
89 template<typename VectorsType, typename CoeffsType, int Side>
90 struct hseq_side_dependent_impl
91 {
92   typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
93   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
94   typedef typename VectorsType::Index Index;
95   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
96   {
97     Index start = k+1+h.m_shift;
98     return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
99   }
100 };
101
102 template<typename VectorsType, typename CoeffsType>
103 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
104 {
105   typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
106   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
107   typedef typename VectorsType::Index Index;
108   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
109   {
110     Index start = k+1+h.m_shift;
111     return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
112   }
113 };
114
115 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
116 {
117   typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
118     ResultScalar;
119   typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
120                  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
121 };
122
123 } // end namespace internal
124
125 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
126   : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
127 {
128     enum {
129       RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
130       ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
131       MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
132       MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
133     };
134     typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
135     typedef typename VectorsType::Index Index;
136
137     typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
138             EssentialVectorType;
139
140   public:
141
142     typedef HouseholderSequence<
143       VectorsType,
144       typename internal::conditional<NumTraits<Scalar>::IsComplex,
145         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
146         CoeffsType>::type,
147       Side
148     > ConjugateReturnType;
149
150     /** \brief Constructor.
151       * \param[in]  v      %Matrix containing the essential parts of the Householder vectors
152       * \param[in]  h      Vector containing the Householder coefficients
153       *
154       * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The
155       * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th
156       * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the
157       * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many
158       * Householder reflections as there are columns.
159       *
160       * \note The %HouseholderSequence object stores \p v and \p h by reference.
161       *
162       * Example: \include HouseholderSequence_HouseholderSequence.cpp
163       * Output: \verbinclude HouseholderSequence_HouseholderSequence.out
164       *
165       * \sa setLength(), setShift()
166       */
167     HouseholderSequence(const VectorsType& v, const CoeffsType& h)
168       : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
169         m_shift(0)
170     {
171     }
172
173     /** \brief Copy constructor. */
174     HouseholderSequence(const HouseholderSequence& other)
175       : m_vectors(other.m_vectors),
176         m_coeffs(other.m_coeffs),
177         m_trans(other.m_trans),
178         m_length(other.m_length),
179         m_shift(other.m_shift)
180     {
181     }
182
183     /** \brief Number of rows of transformation viewed as a matrix.
184       * \returns Number of rows 
185       * \details This equals the dimension of the space that the transformation acts on.
186       */
187     Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
188
189     /** \brief Number of columns of transformation viewed as a matrix.
190       * \returns Number of columns
191       * \details This equals the dimension of the space that the transformation acts on.
192       */
193     Index cols() const { return rows(); }
194
195     /** \brief Essential part of a Householder vector.
196       * \param[in]  k  Index of Householder reflection
197       * \returns    Vector containing non-trivial entries of k-th Householder vector
198       *
199       * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of
200       * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector
201       * \f[ 
202       * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
203       * \f]
204       * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v
205       * passed to the constructor.
206       *
207       * \sa setShift(), shift()
208       */
209     const EssentialVectorType essentialVector(Index k) const
210     {
211       eigen_assert(k >= 0 && k < m_length);
212       return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
213     }
214
215     /** \brief %Transpose of the Householder sequence. */
216     HouseholderSequence transpose() const
217     {
218       return HouseholderSequence(*this).setTrans(!m_trans);
219     }
220
221     /** \brief Complex conjugate of the Householder sequence. */
222     ConjugateReturnType conjugate() const
223     {
224       return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
225              .setTrans(m_trans)
226              .setLength(m_length)
227              .setShift(m_shift);
228     }
229
230     /** \brief Adjoint (conjugate transpose) of the Householder sequence. */
231     ConjugateReturnType adjoint() const
232     {
233       return conjugate().setTrans(!m_trans);
234     }
235
236     /** \brief Inverse of the Householder sequence (equals the adjoint). */
237     ConjugateReturnType inverse() const { return adjoint(); }
238
239     /** \internal */
240     template<typename DestType> void evalTo(DestType& dst) const
241     {
242       Index vecs = m_length;
243       // FIXME find a way to pass this temporary if the user wants to
244       Matrix<Scalar, DestType::RowsAtCompileTime, 1,
245              AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
246       if(    internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value
247           && internal::extract_data(dst) == internal::extract_data(m_vectors))
248       {
249         // in-place
250         dst.diagonal().setOnes();
251         dst.template triangularView<StrictlyUpper>().setZero();
252         for(Index k = vecs-1; k >= 0; --k)
253         {
254           Index cornerSize = rows() - k - m_shift;
255           if(m_trans)
256             dst.bottomRightCorner(cornerSize, cornerSize)
257             .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
258           else
259             dst.bottomRightCorner(cornerSize, cornerSize)
260               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
261
262           // clear the off diagonal vector
263           dst.col(k).tail(rows()-k-1).setZero();
264         }
265         // clear the remaining columns if needed
266         for(Index k = 0; k<cols()-vecs ; ++k)
267           dst.col(k).tail(rows()-k-1).setZero();
268       }
269       else
270       {
271         dst.setIdentity(rows(), rows());
272         for(Index k = vecs-1; k >= 0; --k)
273         {
274           Index cornerSize = rows() - k - m_shift;
275           if(m_trans)
276             dst.bottomRightCorner(cornerSize, cornerSize)
277             .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
278           else
279             dst.bottomRightCorner(cornerSize, cornerSize)
280               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
281         }
282       }
283     }
284
285     /** \internal */
286     template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
287     {
288       Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
289       for(Index k = 0; k < m_length; ++k)
290       {
291         Index actual_k = m_trans ? m_length-k-1 : k;
292         dst.rightCols(rows()-m_shift-actual_k)
293            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
294       }
295     }
296
297     /** \internal */
298     template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
299     {
300       Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
301       for(Index k = 0; k < m_length; ++k)
302       {
303         Index actual_k = m_trans ? k : m_length-k-1;
304         dst.bottomRows(rows()-m_shift-actual_k)
305            .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
306       }
307     }
308
309     /** \brief Computes the product of a Householder sequence with a matrix.
310       * \param[in]  other  %Matrix being multiplied.
311       * \returns    Expression object representing the product.
312       *
313       * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
314       * and \f$ M \f$ is the matrix \p other.
315       */
316     template<typename OtherDerived>
317     typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
318     {
319       typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
320         res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
321       applyThisOnTheLeft(res);
322       return res;
323     }
324
325     template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
326
327     /** \brief Sets the length of the Householder sequence.
328       * \param [in]  length  New value for the length.
329       *
330       * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set
331       * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that
332       * is smaller. After this function is called, the length equals \p length.
333       *
334       * \sa length()
335       */
336     HouseholderSequence& setLength(Index length)
337     {
338       m_length = length;
339       return *this;
340     }
341
342     /** \brief Sets the shift of the Householder sequence.
343       * \param [in]  shift  New value for the shift.
344       *
345       * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th
346       * column of the matrix \p v passed to the constructor corresponds to the i-th Householder
347       * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}}
348       * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th
349       * Householder reflection.
350       *
351       * \sa shift()
352       */
353     HouseholderSequence& setShift(Index shift)
354     {
355       m_shift = shift;
356       return *this;
357     }
358
359     Index length() const { return m_length; }  /**< \brief Returns the length of the Householder sequence. */
360     Index shift() const { return m_shift; }    /**< \brief Returns the shift of the Householder sequence. */
361
362     /* Necessary for .adjoint() and .conjugate() */
363     template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
364
365   protected:
366
367     /** \brief Sets the transpose flag.
368       * \param [in]  trans  New value of the transpose flag.
369       *
370       * By default, the transpose flag is not set. If the transpose flag is set, then this object represents 
371       * \f$ H^T = H_{n-1}^T \ldots H_1^T H_0^T \f$ instead of \f$ H = H_0 H_1 \ldots H_{n-1} \f$.
372       *
373       * \sa trans()
374       */
375     HouseholderSequence& setTrans(bool trans)
376     {
377       m_trans = trans;
378       return *this;
379     }
380
381     bool trans() const { return m_trans; }     /**< \brief Returns the transpose flag. */
382
383     typename VectorsType::Nested m_vectors;
384     typename CoeffsType::Nested m_coeffs;
385     bool m_trans;
386     Index m_length;
387     Index m_shift;
388 };
389
390 /** \brief Computes the product of a matrix with a Householder sequence.
391   * \param[in]  other  %Matrix being multiplied.
392   * \param[in]  h      %HouseholderSequence being multiplied.
393   * \returns    Expression object representing the product.
394   *
395   * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the
396   * Householder sequence represented by \p h.
397   */
398 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
399 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
400 {
401   typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
402     res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
403   h.applyThisOnTheRight(res);
404   return res;
405 }
406
407 /** \ingroup Householder_Module \householder_module
408   * \brief Convenience function for constructing a Householder sequence. 
409   * \returns A HouseholderSequence constructed from the specified arguments.
410   */
411 template<typename VectorsType, typename CoeffsType>
412 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
413 {
414   return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
415 }
416
417 /** \ingroup Householder_Module \householder_module
418   * \brief Convenience function for constructing a Householder sequence. 
419   * \returns A HouseholderSequence constructed from the specified arguments.
420   * \details This function differs from householderSequence() in that the template argument \p OnTheSide of
421   * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft.
422   */
423 template<typename VectorsType, typename CoeffsType>
424 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
425 {
426   return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
427 }
428
429 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H