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

Collapse All | Expand All

(-)a/Eigen/src/QR/HouseholderQR.h (-42 / +48 lines)
Lines 251-306 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename Link Here
251
}
251
}
252
252
253
/** \internal */
253
/** \internal */
254
template<typename MatrixQR, typename HCoeffs>
254
template<typename MatrixQR, typename HCoeffs,
255
void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
255
  typename MatrixQRScalar = typename MatrixQR::Scalar,
256
                                       typename MatrixQR::Index maxBlockSize=32,
256
  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
257
                                       typename MatrixQR::Scalar* tempData = 0)
257
struct householder_qr_inplace_blocked
258
{
258
{
259
  typedef typename MatrixQR::Index Index;
259
  // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
260
  typedef typename MatrixQR::Scalar Scalar;
260
  static void run(MatrixQR& mat, HCoeffs& hCoeffs,
261
  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
261
      typename MatrixQR::Index maxBlockSize=32,
262
262
      typename MatrixQR::Scalar* tempData = 0)
263
  Index rows = mat.rows();
264
  Index cols = mat.cols();
265
  Index size = (std::min)(rows, cols);
266
267
  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
268
  TempType tempVector;
269
  if(tempData==0)
270
  {
263
  {
271
    tempVector.resize(cols);
264
    typedef typename MatrixQR::Index Index;
272
    tempData = tempVector.data();
265
    typedef typename MatrixQR::Scalar Scalar;
273
  }
266
    typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
274
275
  Index blockSize = (std::min)(maxBlockSize,size);
276
277
  Index k = 0;
278
  for (k = 0; k < size; k += blockSize)
279
  {
280
    Index bs = (std::min)(size-k,blockSize);  // actual size of the block
281
    Index tcols = cols - k - bs;            // trailing columns
282
    Index brows = rows-k;                   // rows of the block
283
267
284
    // partition the matrix:
268
    Index rows = mat.rows();
285
    //        A00 | A01 | A02
269
    Index cols = mat.cols();
286
    // mat  = A10 | A11 | A12
270
    Index size = (std::min)(rows, cols);
287
    //        A20 | A21 | A22
288
    // and performs the qr dec of [A11^T A12^T]^T
289
    // and update [A21^T A22^T]^T using level 3 operations.
290
    // Finally, the algorithm continue on A22
291
271
292
    BlockType A11_21 = mat.block(k,k,brows,bs);
272
    typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
293
    Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
273
    TempType tempVector;
274
    if(tempData==0)
275
    {
276
      tempVector.resize(cols);
277
      tempData = tempVector.data();
278
    }
294
279
295
    householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
280
    Index blockSize = (std::min)(maxBlockSize,size);
296
281
297
    if(tcols)
282
    Index k = 0;
283
    for (k = 0; k < size; k += blockSize)
298
    {
284
    {
299
      BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
285
      Index bs = (std::min)(size-k,blockSize);  // actual size of the block
300
      apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
286
      Index tcols = cols - k - bs;            // trailing columns
287
      Index brows = rows-k;                   // rows of the block
288
289
      // partition the matrix:
290
      //        A00 | A01 | A02
291
      // mat  = A10 | A11 | A12
292
      //        A20 | A21 | A22
293
      // and performs the qr dec of [A11^T A12^T]^T
294
      // and update [A21^T A22^T]^T using level 3 operations.
295
      // Finally, the algorithm continue on A22
296
297
      BlockType A11_21 = mat.block(k,k,brows,bs);
298
      Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
299
300
      householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
301
302
      if(tcols)
303
      {
304
        BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
305
        apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
306
      }
301
    }
307
    }
302
  }
308
  }
303
}
309
};
304
310
305
template<typename _MatrixType, typename Rhs>
311
template<typename _MatrixType, typename Rhs>
306
struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
312
struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
Lines 352-358 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& Link Here
352
358
353
  m_temp.resize(cols);
359
  m_temp.resize(cols);
354
360
355
  internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
361
  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
356
362
357
  m_isInitialized = true;
363
  m_isInitialized = true;
358
  return *this;
364
  return *this;
(-)a/Eigen/src/QR/HouseholderQR_MKL.h (-14 / +16 lines)
Lines 34-61 Link Here
34
#ifndef EIGEN_QR_MKL_H
34
#ifndef EIGEN_QR_MKL_H
35
#define EIGEN_QR_MKL_H
35
#define EIGEN_QR_MKL_H
36
36
37
#include "Eigen/src/Core/util/MKL_support.h"
37
#include "../Core/util/MKL_support.h"
38
38
39
namespace Eigen { 
39
namespace Eigen { 
40
40
41
namespace internal {
41
  namespace internal {
42
42
43
/** \internal Specialization for the data types supported by MKL */
43
    /** \internal Specialization for the data types supported by MKL */
44
44
45
#define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \
45
#define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \
46
template<typename MatrixQR, typename HCoeffs> \
46
template<typename MatrixQR, typename HCoeffs> \
47
void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, \
47
struct householder_qr_inplace_blocked<MatrixQR, HCoeffs, EIGTYPE, true> \
48
                                       typename MatrixQR::Index maxBlockSize=32, \
49
                                       EIGTYPE* tempData = 0) \
50
{ \
48
{ \
51
  lapack_int m = mat.rows(); \
49
  static void run(MatrixQR& mat, HCoeffs& hCoeffs, \
52
  lapack_int n = mat.cols(); \
50
      typename MatrixQR::Index = 32, \
53
  lapack_int lda = mat.outerStride(); \
51
      typename MatrixQR::Scalar* = 0) \
54
  lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
52
  { \
55
  LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \
53
    lapack_int m = (lapack_int) mat.rows(); \
56
  hCoeffs.adjointInPlace(); \
54
    lapack_int n = (lapack_int) mat.cols(); \
57
\
55
    lapack_int lda = (lapack_int) mat.outerStride(); \
58
}
56
    lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
57
    LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \
58
    hCoeffs.adjointInPlace(); \
59
  } \
60
};
59
61
60
EIGEN_MKL_QR_NOPIV(double, double, d)
62
EIGEN_MKL_QR_NOPIV(double, double, d)
61
EIGEN_MKL_QR_NOPIV(float, float, s)
63
EIGEN_MKL_QR_NOPIV(float, float, s)

Return to bug 704