Jaesung Oh / Eigen

Eigne Matrix Class Library

Dependents:   MPC_current_control HydraulicControlBoard_SW AHRS Test_ekf ... more

src/Householder/HouseholderSequence.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

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