diff options
author | karturov <karturov@KARTUROV-MOBL.ccr.corp.intel.com> | 2011-12-05 14:52:21 +0700 |
---|---|---|
committer | karturov <karturov@KARTUROV-MOBL.ccr.corp.intel.com> | 2011-12-05 14:52:21 +0700 |
commit | 015c331252a3b99c187b5607572f1cec531a4d1e (patch) | |
tree | e30a3f64a950edd21ae89f667cb2f859b6479c02 /Eigen/src/Cholesky/LLT_MKL.h | |
parent | e270a5656aaafb055702a51a63541e05eabd8936 (diff) |
Intel(R) MKL support added.
* * *
License disclaimer changed to BSD license for MKL_support.h
* * *
Pardiso support fixed, test added.
blas/lapack tests fixed: Scalar parameter was added in Cholesky, product_matrix_vector_triangular remaned to triangular_matrix_vector_product.
* * *
PARDISO test was added physically.
Diffstat (limited to 'Eigen/src/Cholesky/LLT_MKL.h')
-rw-r--r-- | Eigen/src/Cholesky/LLT_MKL.h | 123 |
1 files changed, 123 insertions, 0 deletions
diff --git a/Eigen/src/Cholesky/LLT_MKL.h b/Eigen/src/Cholesky/LLT_MKL.h new file mode 100644 index 000000000..c5ed77d33 --- /dev/null +++ b/Eigen/src/Cholesky/LLT_MKL.h @@ -0,0 +1,123 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to Intel(R) MKL + * LLt decomposition based on LAPACKE_?potrf function. + ******************************************************************************** +*/ + +#ifndef EIGEN_LLT_MKL_H +#define EIGEN_LLT_MKL_H + +#include "Eigen/src/Core/util/MKL_support.h" +#include <iostream> + +namespace internal { + +template<typename Scalar> struct mkl_llt; + +#define EIGEN_MKL_LLT(EIGTYPE, MKLTYPE, MKLPREFIX) \ +template<> struct mkl_llt<EIGTYPE> \ +{ \ + template<typename MatrixType> \ + static inline typename MatrixType::Index potrf(MatrixType& m, char uplo) \ + { \ + lapack_int matrix_order; \ + lapack_int size, lda, info, StorageOrder; \ + EIGTYPE* a; \ + eigen_assert(m.rows()==m.cols()); \ +/* Set up parameters for ?potrf */ \ + size = m.rows(); \ + StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \ + matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ + a = &(m.coeffRef(0,0)); \ + lda = m.outerStride(); \ +\ + info = LAPACKE_##MKLPREFIX##potrf( matrix_order, uplo, size, (MKLTYPE*)a, lda ); \ + info = (info==0) ? Success : NumericalIssue; \ + return info; \ + } \ +}; \ +template<> struct llt_inplace<EIGTYPE, Lower> \ +{ \ + template<typename MatrixType> \ + static typename MatrixType::Index blocked(MatrixType& m) \ + { \ + return mkl_llt<EIGTYPE>::potrf(m, 'L'); \ + } \ + template<typename MatrixType, typename VectorType> \ + static void rankUpdate(MatrixType& mat, const VectorType& vec) \ + { \ + typedef typename MatrixType::ColXpr ColXpr; \ + typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; \ + typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef Matrix<Scalar,Dynamic,1> TempVectorType; \ + typedef typename TempVectorType::SegmentReturnType TempVecSegment; \ +\ + int n = mat.cols(); \ + eigen_assert(mat.rows()==n && vec.size()==n); \ + TempVectorType temp(vec); \ +\ + for(int i=0; i<n; ++i) \ + { \ + JacobiRotation<Scalar> g; \ + g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); \ +\ + int rs = n-i-1; \ + if(rs>0) \ + { \ + ColXprSegment x(mat.col(i).tail(rs)); \ + TempVecSegment y(temp.tail(rs)); \ + apply_rotation_in_the_plane(x, y, g); \ + } \ + } \ + } \ +}; \ +template<> struct llt_inplace<EIGTYPE, Upper> \ +{ \ + template<typename MatrixType> \ + static typename MatrixType::Index blocked(MatrixType& m) \ + { \ + return mkl_llt<EIGTYPE>::potrf(m, 'U'); \ + } \ + template<typename MatrixType, typename VectorType> \ + static void rankUpdate(MatrixType& mat, const VectorType& vec) \ + { \ + Transpose<MatrixType> matt(mat); \ + return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate()); \ + } \ +}; + +EIGEN_MKL_LLT(double, double, d) +EIGEN_MKL_LLT(float, float, s) +EIGEN_MKL_LLT(dcomplex, MKL_Complex16, z) +EIGEN_MKL_LLT(scomplex, MKL_Complex8, c) + +} + +#endif // EIGEN_LLT_MKL_H |