aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/Cholesky/LLT_MKL.h
diff options
context:
space:
mode:
authorGravatar karturov <karturov@KARTUROV-MOBL.ccr.corp.intel.com>2011-12-05 14:52:21 +0700
committerGravatar karturov <karturov@KARTUROV-MOBL.ccr.corp.intel.com>2011-12-05 14:52:21 +0700
commit015c331252a3b99c187b5607572f1cec531a4d1e (patch)
treee30a3f64a950edd21ae89f667cb2f859b6479c02 /Eigen/src/Cholesky/LLT_MKL.h
parente270a5656aaafb055702a51a63541e05eabd8936 (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.h123
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