Applied patch from Richard JW Roberts, resolving bug #704

(grafted from devel branch)
This commit is contained in:
Gael Guennebaud 2015-06-15 22:02:57 +02:00
parent 1c6b224fb3
commit 595c00157c
2 changed files with 66 additions and 58 deletions

View File

@ -257,8 +257,13 @@ void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename
} }
/** \internal */ /** \internal */
template<typename MatrixQR, typename HCoeffs> template<typename MatrixQR, typename HCoeffs,
void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQRScalar = typename MatrixQR::Scalar,
bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
struct householder_qr_inplace_blocked
{
// This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
static void run(MatrixQR& mat, HCoeffs& hCoeffs,
typename MatrixQR::Index maxBlockSize=32, typename MatrixQR::Index maxBlockSize=32,
typename MatrixQR::Scalar* tempData = 0) typename MatrixQR::Scalar* tempData = 0)
{ {
@ -307,6 +312,7 @@ void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
} }
} }
} }
};
template<typename _MatrixType, typename Rhs> template<typename _MatrixType, typename Rhs>
struct solve_retval<HouseholderQR<_MatrixType>, Rhs> struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
@ -360,7 +366,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
m_temp.resize(cols); m_temp.resize(cols);
internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data()); internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;

View File

@ -34,7 +34,7 @@
#ifndef EIGEN_QR_MKL_H #ifndef EIGEN_QR_MKL_H
#define EIGEN_QR_MKL_H #define EIGEN_QR_MKL_H
#include "Eigen/src/Core/util/MKL_support.h" #include "../Core/util/MKL_support.h"
namespace Eigen { namespace Eigen {
@ -44,18 +44,20 @@ namespace internal {
#define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \ #define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \
template<typename MatrixQR, typename HCoeffs> \ template<typename MatrixQR, typename HCoeffs> \
void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, \ struct householder_qr_inplace_blocked<MatrixQR, HCoeffs, EIGTYPE, true> \
typename MatrixQR::Index maxBlockSize=32, \
EIGTYPE* tempData = 0) \
{ \ { \
lapack_int m = mat.rows(); \ static void run(MatrixQR& mat, HCoeffs& hCoeffs, \
lapack_int n = mat.cols(); \ typename MatrixQR::Index = 32, \
lapack_int lda = mat.outerStride(); \ typename MatrixQR::Scalar* = 0) \
{ \
lapack_int m = (lapack_int) mat.rows(); \
lapack_int n = (lapack_int) mat.cols(); \
lapack_int lda = (lapack_int) mat.outerStride(); \
lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \ LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \
hCoeffs.adjointInPlace(); \ hCoeffs.adjointInPlace(); \
\ } \
} };
EIGEN_MKL_QR_NOPIV(double, double, d) EIGEN_MKL_QR_NOPIV(double, double, d)
EIGEN_MKL_QR_NOPIV(float, float, s) EIGEN_MKL_QR_NOPIV(float, float, s)