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

Collapse All | Expand All

(-)a/Eigen/src/Householder/HouseholderSequence.h (-13 / +36 lines)
Lines 232-313 template<typename VectorsType, typename Link Here
232
    {
232
    {
233
      return conjugate().setTrans(!m_trans);
233
      return conjugate().setTrans(!m_trans);
234
    }
234
    }
235
235
236
    /** \brief Inverse of the Householder sequence (equals the adjoint). */
236
    /** \brief Inverse of the Householder sequence (equals the adjoint). */
237
    ConjugateReturnType inverse() const { return adjoint(); }
237
    ConjugateReturnType inverse() const { return adjoint(); }
238
238
239
    /** \internal */
239
    /** \internal */
240
    template<typename DestType> void evalTo(DestType& dst) const
240
    template<typename DestType> inline void evalTo(DestType& dst) const
241
    {
241
    {
242
      Matrix<Scalar, DestType::RowsAtCompileTime, 1,
243
             AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
244
      evalTo(dst, workspace);
245
    }
246
247
    /** \internal */
248
    template<typename Dest, typename Workspace>
249
    void evalTo(Dest& dst, MatrixBase<Workspace>& workspace) const
250
    {
251
      workspace.derived().resize(rows());
242
      Index vecs = m_length;
252
      Index vecs = m_length;
243
      // FIXME find a way to pass this temporary if the user wants to
253
      if(    internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
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))
254
          && internal::extract_data(dst) == internal::extract_data(m_vectors))
248
      {
255
      {
249
        // in-place
256
        // in-place
250
        dst.diagonal().setOnes();
257
        dst.diagonal().setOnes();
251
        dst.template triangularView<StrictlyUpper>().setZero();
258
        dst.template triangularView<StrictlyUpper>().setZero();
252
        for(Index k = vecs-1; k >= 0; --k)
259
        for(Index k = vecs-1; k >= 0; --k)
253
        {
260
        {
254
          Index cornerSize = rows() - k - m_shift;
261
          Index cornerSize = rows() - k - m_shift;
255
          if(m_trans)
262
          if(m_trans)
256
            dst.bottomRightCorner(cornerSize, cornerSize)
263
            dst.bottomRightCorner(cornerSize, cornerSize)
257
            .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
264
               .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
258
          else
265
          else
259
            dst.bottomRightCorner(cornerSize, cornerSize)
266
            dst.bottomRightCorner(cornerSize, cornerSize)
260
              .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
267
               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
261
268
262
          // clear the off diagonal vector
269
          // clear the off diagonal vector
263
          dst.col(k).tail(rows()-k-1).setZero();
270
          dst.col(k).tail(rows()-k-1).setZero();
264
        }
271
        }
265
        // clear the remaining columns if needed
272
        // clear the remaining columns if needed
266
        for(Index k = 0; k<cols()-vecs ; ++k)
273
        for(Index k = 0; k<cols()-vecs ; ++k)
267
          dst.col(k).tail(rows()-k-1).setZero();
274
          dst.col(k).tail(rows()-k-1).setZero();
268
      }
275
      }
269
      else
276
      else
270
      {
277
      {
271
        dst.setIdentity(rows(), rows());
278
        dst.setIdentity(rows(), rows());
272
        for(Index k = vecs-1; k >= 0; --k)
279
        for(Index k = vecs-1; k >= 0; --k)
273
        {
280
        {
274
          Index cornerSize = rows() - k - m_shift;
281
          Index cornerSize = rows() - k - m_shift;
275
          if(m_trans)
282
          if(m_trans)
276
            dst.bottomRightCorner(cornerSize, cornerSize)
283
            dst.bottomRightCorner(cornerSize, cornerSize)
277
            .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
284
               .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
278
          else
285
          else
279
            dst.bottomRightCorner(cornerSize, cornerSize)
286
            dst.bottomRightCorner(cornerSize, cornerSize)
280
              .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
287
               .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
281
        }
288
        }
282
      }
289
      }
283
    }
290
    }
284
291
285
    /** \internal */
292
    /** \internal */
286
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
293
    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
287
    {
294
    {
288
      Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
295
      Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
296
      applyThisOnTheRight(dst, workspace);
297
    }
298
299
    /** \internal */
300
    template<typename Dest, typename Workspace>
301
    inline void applyThisOnTheRight(Dest& dst, MatrixBase<Workspace>& workspace) const
302
    {
303
      workspace.derived().resize(dst.rows());
289
      for(Index k = 0; k < m_length; ++k)
304
      for(Index k = 0; k < m_length; ++k)
290
      {
305
      {
291
        Index actual_k = m_trans ? m_length-k-1 : k;
306
        Index actual_k = m_trans ? m_length-k-1 : k;
292
        dst.rightCols(rows()-m_shift-actual_k)
307
        dst.rightCols(rows()-m_shift-actual_k)
293
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
308
           .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &workspace.coeffRef(0));
294
      }
309
      }
295
    }
310
    }
296
311
297
    /** \internal */
312
    /** \internal */
298
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
313
    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
299
    {
314
    {
300
      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
315
      Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
316
      applyThisOnTheLeft(dst, workspace);
317
    }
318
319
    /** \internal */
320
    template<typename Dest, typename Workspace>
321
    inline void applyThisOnTheLeft(Dest& dst, MatrixBase<Workspace>& workspace) const
322
    {
323
      workspace.derived().resize(dst.cols());
301
      for(Index k = 0; k < m_length; ++k)
324
      for(Index k = 0; k < m_length; ++k)
302
      {
325
      {
303
        Index actual_k = m_trans ? k : m_length-k-1;
326
        Index actual_k = m_trans ? k : m_length-k-1;
304
        dst.bottomRows(rows()-m_shift-actual_k)
327
        dst.bottomRows(rows()-m_shift-actual_k)
305
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
328
           .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &workspace.coeffRef(0));
306
      }
329
      }
307
    }
330
    }
308
331
309
    /** \brief Computes the product of a Householder sequence with a matrix.
332
    /** \brief Computes the product of a Householder sequence with a matrix.
310
      * \param[in]  other  %Matrix being multiplied.
333
      * \param[in]  other  %Matrix being multiplied.
311
      * \returns    Expression object representing the product.
334
      * \returns    Expression object representing the product.
312
      *
335
      *
313
      * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
336
      * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this

Return to bug 206