1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*
2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Copyright (c) 2011, Intel Corporation. All rights reserved.
3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Redistribution and use in source and binary forms, with or without modification,
5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath are permitted provided that the following conditions are met:
6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Redistributions of source code must retain the above copyright notice, this
8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath   list of conditions and the following disclaimer.
9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Redistributions in binary form must reproduce the above copyright notice,
10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath   this list of conditions and the following disclaimer in the documentation
11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath   and/or other materials provided with the distribution.
12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Neither the name of Intel Corporation nor the names of its contributors may
13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath   be used to endorse or promote products derived from this software without
14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath   specific prior written permission.
15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ********************************************************************************
28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *   Content : Eigen bindings to Intel(R) MKL
29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *     LLt decomposition based on LAPACKE_?potrf function.
30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ********************************************************************************
31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath*/
32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_LLT_MKL_H
34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_LLT_MKL_H
35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "Eigen/src/Core/util/MKL_support.h"
37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <iostream>
38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen {
40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal {
42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> struct mkl_llt;
44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_MKL_LLT(EIGTYPE, MKLTYPE, MKLPREFIX) \
46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<> struct mkl_llt<EIGTYPE> \
47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ \
48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  template<typename MatrixType> \
49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static inline typename MatrixType::Index potrf(MatrixType& m, char uplo) \
50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  { \
51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    lapack_int matrix_order; \
52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    lapack_int size, lda, info, StorageOrder; \
53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    EIGTYPE* a; \
54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    eigen_assert(m.rows()==m.cols()); \
55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    /* Set up parameters for ?potrf */ \
56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    size = m.rows(); \
57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \
58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    a = &(m.coeffRef(0,0)); \
60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    lda = m.outerStride(); \
61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath\
62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    info = LAPACKE_##MKLPREFIX##potrf( matrix_order, uplo, size, (MKLTYPE*)a, lda ); \
63a829215e078ace896f52702caa0c27608f40e3b0Miao Wang    info = (info==0) ? -1 : info>0 ? info-1 : size; \
64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    return info; \
65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  } \
66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; \
67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<> struct llt_inplace<EIGTYPE, Lower> \
68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ \
69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  template<typename MatrixType> \
70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static typename MatrixType::Index blocked(MatrixType& m) \
71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  { \
72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    return mkl_llt<EIGTYPE>::potrf(m, 'L'); \
73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  } \
74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  template<typename MatrixType, typename VectorType> \
75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  { return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); } \
77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; \
78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<> struct llt_inplace<EIGTYPE, Upper> \
79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ \
80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  template<typename MatrixType> \
81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static typename MatrixType::Index blocked(MatrixType& m) \
82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  { \
83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    return mkl_llt<EIGTYPE>::potrf(m, 'U'); \
84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  } \
85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  template<typename MatrixType, typename VectorType> \
86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  { \
88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    Transpose<MatrixType> matt(mat); \
89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate(), sigma); \
90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  } \
91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath};
92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_MKL_LLT(double, double, d)
94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_MKL_LLT(float, float, s)
95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_MKL_LLT(dcomplex, MKL_Complex16, z)
96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_MKL_LLT(scomplex, MKL_Complex8, c)
97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal
99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen
101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_LLT_MKL_H
103