This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 50 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/Eigenvalues/HessenbergDecomposition.h (-1 / +3 lines)
Lines 240-256 template<typename _MatrixType> class Hes Link Here
240
      * HouseholderSequence. You can either apply it directly to a matrix or
240
      * HouseholderSequence. You can either apply it directly to a matrix or
241
      * you can convert it to a matrix of type #MatrixType.
241
      * you can convert it to a matrix of type #MatrixType.
242
      *
242
      *
243
      * \sa matrixH() for an example, class HouseholderSequence
243
      * \sa matrixH() for an example, class HouseholderSequence
244
      */
244
      */
245
    HouseholderSequenceType matrixQ() const
245
    HouseholderSequenceType matrixQ() const
246
    {
246
    {
247
      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
247
      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
248
      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
248
      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
249
             .setLength(m_matrix.rows() - 1)
250
             .setShift(1);
249
    }
251
    }
250
252
251
    /** \brief Constructs the Hessenberg matrix H in the decomposition
253
    /** \brief Constructs the Hessenberg matrix H in the decomposition
252
      *
254
      *
253
      * \returns expression object representing the matrix H
255
      * \returns expression object representing the matrix H
254
      *
256
      *
255
      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
257
      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
256
      * or the member function compute(const MatrixType&) has been called
258
      * or the member function compute(const MatrixType&) has been called
(-)a/Eigen/src/Eigenvalues/Tridiagonalization.h (-2 / +6 lines)
Lines 246-262 template<typename _MatrixType> class Tri Link Here
246
      * you can convert it to a matrix of type #MatrixType.
246
      * you can convert it to a matrix of type #MatrixType.
247
      *
247
      *
248
      * \sa Tridiagonalization(const MatrixType&) for an example,
248
      * \sa Tridiagonalization(const MatrixType&) for an example,
249
      *     matrixT(), class HouseholderSequence
249
      *     matrixT(), class HouseholderSequence
250
      */
250
      */
251
    HouseholderSequenceType matrixQ() const
251
    HouseholderSequenceType matrixQ() const
252
    {
252
    {
253
      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
253
      eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
254
      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
254
      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
255
             .setLength(m_matrix.rows() - 1)
256
             .setShift(1);
255
    }
257
    }
256
258
257
    /** \brief Returns an expression of the tridiagonal matrix T in the decomposition
259
    /** \brief Returns an expression of the tridiagonal matrix T in the decomposition
258
      *
260
      *
259
      * \returns expression object representing the matrix T
261
      * \returns expression object representing the matrix T
260
      *
262
      *
261
      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
263
      * \pre Either the constructor Tridiagonalization(const MatrixType&) or
262
      * the member function compute(const MatrixType&) has been called before
264
      * the member function compute(const MatrixType&) has been called before
Lines 454-470 struct tridiagonalization_inplace_select Link Here
454
  template<typename DiagonalType, typename SubDiagonalType>
456
  template<typename DiagonalType, typename SubDiagonalType>
455
  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
457
  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
456
  {
458
  {
457
    CoeffVectorType hCoeffs(mat.cols()-1);
459
    CoeffVectorType hCoeffs(mat.cols()-1);
458
    tridiagonalization_inplace(mat,hCoeffs);
460
    tridiagonalization_inplace(mat,hCoeffs);
459
    diag = mat.diagonal().real();
461
    diag = mat.diagonal().real();
460
    subdiag = mat.template diagonal<-1>().real();
462
    subdiag = mat.template diagonal<-1>().real();
461
    if(extractQ)
463
    if(extractQ)
462
      mat = HouseholderSequenceType(mat, hCoeffs.conjugate(), false, mat.rows() - 1, 1);
464
      mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
465
            .setLength(mat.rows() - 1)
466
            .setShift(1);
463
  }
467
  }
464
};
468
};
465
469
466
/** \internal
470
/** \internal
467
  * Specialization for 3x3 real matrices.
471
  * Specialization for 3x3 real matrices.
468
  * Especially useful for plane fitting.
472
  * Especially useful for plane fitting.
469
  */
473
  */
470
template<typename MatrixType>
474
template<typename MatrixType>
(-)a/Eigen/src/Householder/HouseholderSequence.h (-34 / +53 lines)
Lines 124-174 template<typename VectorsType, typename Link Here
124
    typedef HouseholderSequence<
124
    typedef HouseholderSequence<
125
      VectorsType,
125
      VectorsType,
126
      typename internal::conditional<NumTraits<Scalar>::IsComplex,
126
      typename internal::conditional<NumTraits<Scalar>::IsComplex,
127
        typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
127
        typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
128
        CoeffsType>::type,
128
        CoeffsType>::type,
129
      Side
129
      Side
130
    > ConjugateReturnType;
130
    > ConjugateReturnType;
131
131
132
    HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
132
    HouseholderSequence(const VectorsType& v, const CoeffsType& h)
133
      : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()),
133
      : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
134
        m_shift(0)
134
        m_shift(0)
135
    {
135
    {
136
    }
136
    }
137
137
138
    HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, Index actualVectors, Index shift)
138
    HouseholderSequence(const HouseholderSequence& other)
139
      : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors), m_shift(shift)
139
      : m_vectors(other.m_vectors),
140
        m_coeffs(other.m_coeffs),
141
        m_trans(other.m_trans),
142
        m_length(other.m_length),
143
        m_shift(other.m_shift)
140
    {
144
    {
141
    }
145
    }
142
146
143
    Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
147
    Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
144
    Index cols() const { return rows(); }
148
    Index cols() const { return rows(); }
145
149
146
    const EssentialVectorType essentialVector(Index k) const
150
    const EssentialVectorType essentialVector(Index k) const
147
    {
151
    {
148
      eigen_assert(k >= 0 && k < m_actualVectors);
152
      eigen_assert(k >= 0 && k < m_length);
149
      return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
153
      return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
150
    }
154
    }
151
155
152
    HouseholderSequence transpose() const
156
    HouseholderSequence transpose() const
153
    { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors, m_shift); }
157
    {
158
      return HouseholderSequence(*this).setTrans(!m_trans);
159
    }
154
160
155
    ConjugateReturnType conjugate() const
161
    ConjugateReturnType conjugate() const
156
    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors, m_shift); }
162
    {
163
      return ConjugateReturnType(m_vectors, m_coeffs.conjugate())
164
             .setTrans(m_trans)
165
             .setLength(m_length)
166
             .setShift(m_shift);
167
    }
157
168
158
    ConjugateReturnType adjoint() const
169
    ConjugateReturnType adjoint() const
159
    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors, m_shift); }
170
    {
171
      return conjugate().setTrans(!m_trans);
172
    }
160
173
161
    ConjugateReturnType inverse() const { return adjoint(); }
174
    ConjugateReturnType inverse() const { return adjoint(); }
162
175
163
    /** \internal */
176
    /** \internal */
164
    template<typename DestType> void evalTo(DestType& dst) const
177
    template<typename DestType> void evalTo(DestType& dst) const
165
    {
178
    {
166
      Index vecs = m_actualVectors;
179
      Index vecs = m_length;
167
      // FIXME find a way to pass this temporary if the user wants to
180
      // FIXME find a way to pass this temporary if the user wants to
168
      Matrix<Scalar, DestType::RowsAtCompileTime, 1,
181
      Matrix<Scalar, DestType::RowsAtCompileTime, 1,
169
             AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
182
             AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> temp(rows());
170
      if(    internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value
183
      if(    internal::is_same<typename internal::remove_all<VectorsType>::type,DestType>::value
171
          && internal::extract_data(dst) == internal::extract_data(m_vectors))
184
          && internal::extract_data(dst) == internal::extract_data(m_vectors))
172
      {
185
      {
173
        // in-place
186
        // in-place
174
        dst.diagonal().setOnes();
187
        dst.diagonal().setOnes();
Lines 205-235 template<typename VectorsType, typename Link Here
205
        }
218
        }
206
      }
219
      }
207
    }
220
    }
208
221
209
    /** \internal */
222
    /** \internal */
210
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
223
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
211
    {
224
    {
212
      Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
225
      Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
213
      for(Index k = 0; k < m_actualVectors; ++k)
226
      for(Index k = 0; k < m_length; ++k)
214
      {
227
      {
215
        Index actual_k = m_trans ? m_actualVectors-k-1 : k;
228
        Index actual_k = m_trans ? m_length-k-1 : k;
216
        dst.rightCols(rows()-m_shift-actual_k)
229
        dst.rightCols(rows()-m_shift-actual_k)
217
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
230
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
218
      }
231
      }
219
    }
232
    }
220
233
221
    /** \internal */
234
    /** \internal */
222
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
235
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
223
    {
236
    {
224
      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
237
      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
225
      for(Index k = 0; k < m_actualVectors; ++k)
238
      for(Index k = 0; k < m_length; ++k)
226
      {
239
      {
227
        Index actual_k = m_trans ? k : m_actualVectors-k-1;
240
        Index actual_k = m_trans ? k : m_length-k-1;
228
        dst.bottomRows(rows()-m_shift-actual_k)
241
        dst.bottomRows(rows()-m_shift-actual_k)
229
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
242
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
230
      }
243
      }
231
    }
244
    }
232
245
233
    template<typename OtherDerived>
246
    template<typename OtherDerived>
234
    typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
247
    typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
235
    {
248
    {
Lines 245-289 template<typename VectorsType, typename Link Here
245
      typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
258
      typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
246
        res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::ResultScalar>());
259
        res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::ResultScalar>());
247
      h.applyThisOnTheRight(res);
260
      h.applyThisOnTheRight(res);
248
      return res;
261
      return res;
249
    }
262
    }
250
263
251
    template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
264
    template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
252
265
266
    HouseholderSequence& setTrans(bool trans)
267
    {
268
      m_trans = trans;
269
      return *this;
270
    }
271
272
    HouseholderSequence& setLength(Index length)
273
    {
274
      m_length = length;
275
      return *this;
276
    }
277
278
    HouseholderSequence& setShift(Index shift)
279
    {
280
      m_shift = shift;
281
      return *this;
282
    }
283
284
    bool trans() const { return m_trans; }
285
    Index length() const { return m_length; }
286
    Index shift() const { return m_shift; }
287
253
  protected:
288
  protected:
254
    typename VectorsType::Nested m_vectors;
289
    typename VectorsType::Nested m_vectors;
255
    typename CoeffsType::Nested m_coeffs;
290
    typename CoeffsType::Nested m_coeffs;
256
    bool m_trans;
291
    bool m_trans;
257
    Index m_actualVectors;
292
    Index m_length;
258
    Index m_shift;
293
    Index m_shift;
259
};
294
};
260
295
261
template<typename VectorsType, typename CoeffsType>
296
template<typename VectorsType, typename CoeffsType>
262
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
297
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
263
{
298
{
264
  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans);
299
  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
265
}
300
}
266
301
267
template<typename VectorsType, typename CoeffsType>
302
template<typename VectorsType, typename CoeffsType>
268
HouseholderSequence<VectorsType,CoeffsType> householderSequence
303
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
269
   (const VectorsType& v, const CoeffsType& h,
270
    bool trans, typename VectorsType::Index actualVectors, typename VectorsType::Index shift)
271
{
304
{
272
  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans, actualVectors, shift);
305
  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
273
}
274
275
template<typename VectorsType, typename CoeffsType>
276
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
277
{
278
  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans);
279
}
280
281
template<typename VectorsType, typename CoeffsType>
282
HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence
283
  (const VectorsType& v, const CoeffsType& h, bool trans,
284
   typename VectorsType::Index actualVectors, typename VectorsType::Index shift)
285
{
286
  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans, actualVectors, shift);
287
}
306
}
288
307
289
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
308
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
(-)a/Eigen/src/QR/ColPivHouseholderQR.h (-8 / +5 lines)
Lines 478-500 struct solve_retval<ColPivHouseholderQR< Link Here
478
    {
478
    {
479
      dst.setZero();
479
      dst.setZero();
480
      return;
480
      return;
481
    }
481
    }
482
482
483
    typename Rhs::PlainObject c(rhs());
483
    typename Rhs::PlainObject c(rhs());
484
484
485
    // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
485
    // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
486
    c.applyOnTheLeft(householderSequence(
486
    c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
487
      dec().matrixQR(),
487
                     .setTrans(true)
488
      dec().hCoeffs(),
488
                     .setLength(dec().nonzeroPivots())
489
      true,
489
      );
490
      dec().nonzeroPivots(),
491
      0
492
    ));
493
490
494
    dec().matrixQR()
491
    dec().matrixQR()
495
       .topLeftCorner(nonzero_pivots, nonzero_pivots)
492
       .topLeftCorner(nonzero_pivots, nonzero_pivots)
496
       .template triangularView<Upper>()
493
       .template triangularView<Upper>()
497
       .solveInPlace(c.topRows(nonzero_pivots));
494
       .solveInPlace(c.topRows(nonzero_pivots));
498
495
499
496
500
    typename Rhs::PlainObject d(c);
497
    typename Rhs::PlainObject d(c);
Lines 512-528 struct solve_retval<ColPivHouseholderQR< Link Here
512
} // end namespace internal
509
} // end namespace internal
513
510
514
/** \returns the matrix Q as a sequence of householder transformations */
511
/** \returns the matrix Q as a sequence of householder transformations */
515
template<typename MatrixType>
512
template<typename MatrixType>
516
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
513
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
517
  ::householderQ() const
514
  ::householderQ() const
518
{
515
{
519
  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
516
  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
520
  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots, 0);
517
  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
521
}
518
}
522
519
523
/** \return the column-pivoting Householder QR decomposition of \c *this.
520
/** \return the column-pivoting Householder QR decomposition of \c *this.
524
  *
521
  *
525
  * \sa class ColPivHouseholderQR
522
  * \sa class ColPivHouseholderQR
526
  */
523
  */
527
template<typename Derived>
524
template<typename Derived>
528
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
525
const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
(-)a/Eigen/src/SVD/UpperBidiagonalization.h (-2 / +3 lines)
Lines 82-99 template<typename _MatrixType> class Upp Link Here
82
    {
82
    {
83
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
83
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
84
      return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
84
      return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
85
    }
85
    }
86
86
87
    const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
87
    const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
88
    {
88
    {
89
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
89
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
90
      return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>(),
90
      return HouseholderVSequenceType(m_householder, m_householder.const_derived().template diagonal<1>())
91
                                      false, m_householder.cols()-1, 1);
91
             .setLength(m_householder.cols()-1)
92
             .setShift(1);
92
    }
93
    }
93
    
94
    
94
  protected:
95
  protected:
95
    MatrixType m_householder;
96
    MatrixType m_householder;
96
    BidiagonalType m_bidiagonal;
97
    BidiagonalType m_bidiagonal;
97
    bool m_isInitialized;
98
    bool m_isInitialized;
98
};
99
};
99
100
(-)a/test/householder.cpp (-2 / +8 lines)
Lines 97-123 template<typename MatrixType> void house Link Here
97
  Index shift = internal::random<Index>(0, std::max<Index>(rows-2,0));
97
  Index shift = internal::random<Index>(0, std::max<Index>(rows-2,0));
98
  Index brows = rows - shift;
98
  Index brows = rows - shift;
99
  m1.setRandom(rows, cols);
99
  m1.setRandom(rows, cols);
100
  HBlockMatrixType hbm = m1.block(shift,0,brows,cols);
100
  HBlockMatrixType hbm = m1.block(shift,0,brows,cols);
101
  HouseholderQR<HBlockMatrixType> qr(hbm);
101
  HouseholderQR<HBlockMatrixType> qr(hbm);
102
  m2 = m1;
102
  m2 = m1;
103
  m2.block(shift,0,brows,cols) = qr.matrixQR();
103
  m2.block(shift,0,brows,cols) = qr.matrixQR();
104
  HCoeffsVectorType hc = qr.hCoeffs().conjugate();
104
  HCoeffsVectorType hc = qr.hCoeffs().conjugate();
105
  HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift);
105
  HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc);
106
  hseq.setLength(hc.size()).setShift(shift);
107
  VERIFY(hseq.trans() == false);
108
  VERIFY(hseq.length() == hc.size());
109
  VERIFY(hseq.shift() == shift);
110
106
  MatrixType m5 = m2;
111
  MatrixType m5 = m2;
107
  m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
112
  m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
108
  VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
113
  VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
109
  m3 = hseq;
114
  m3 = hseq;
110
  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying
115
  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying
111
116
112
  // test householder sequence on the right with a shift
117
  // test householder sequence on the right with a shift
113
118
114
  TMatrixType tm2 = m2.transpose();
119
  TMatrixType tm2 = m2.transpose();
115
  HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc, false, hc.size(), shift);
120
  HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc);
121
  rhseq.setLength(hc.size()).setShift(shift);
116
  VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
122
  VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
117
  m3 = rhseq;
123
  m3 = rhseq;
118
  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying
124
  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying
119
}
125
}
120
126
121
void test_householder()
127
void test_householder()
122
{
128
{
123
  for(int i = 0; i < g_repeat; i++) {
129
  for(int i = 0; i < g_repeat; i++) {

Return to bug 50