diff options
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 20 | ||||
-rw-r--r-- | Eigen/src/Cholesky/LLT_MKL.h | 123 |
2 files changed, 133 insertions, 10 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index f061cbf42..aef7ae263 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -196,15 +196,15 @@ template<typename _MatrixType, int _UpLo> class LLT namespace internal { -template<int UpLo> struct llt_inplace; +template<typename Scalar, int UpLo> struct llt_inplace; -template<> struct llt_inplace<Lower> +template<typename Scalar> struct llt_inplace<Scalar, Lower> { template<typename MatrixType> static typename MatrixType::Index unblocked(MatrixType& mat) { typedef typename MatrixType::Index Index; - typedef typename MatrixType::Scalar Scalar; +// typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; eigen_assert(mat.rows()==mat.cols()); @@ -291,25 +291,25 @@ template<> struct llt_inplace<Lower> } }; -template<> struct llt_inplace<Upper> +template<typename Scalar> struct llt_inplace<Scalar, Upper> { template<typename MatrixType> static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return llt_inplace<Lower>::unblocked(matt); + return llt_inplace<Scalar, Lower>::unblocked(matt); } template<typename MatrixType> static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); - return llt_inplace<Lower>::blocked(matt); + return llt_inplace<Scalar, Lower>::blocked(matt); } template<typename MatrixType, typename VectorType> static void rankUpdate(MatrixType& mat, const VectorType& vec) { Transpose<MatrixType> matt(mat); - return llt_inplace<Lower>::rankUpdate(matt, vec.conjugate()); + return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate()); } }; @@ -320,7 +320,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Lower>::blocked(m)==-1; } + { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; } }; template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> @@ -330,7 +330,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Upper>::blocked(m)==-1; } + { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; } }; } // end namespace internal @@ -368,7 +368,7 @@ template<typename VectorType> void LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); - internal::llt_inplace<UpLo>::rankUpdate(m_matrix,v); + internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v); } namespace internal { 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 |