$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 /* 00002 Copyright (c) 2011, Intel Corporation. All rights reserved. 00003 00004 Redistribution and use in source and binary forms, with or without modification, 00005 are permitted provided that the following conditions are met: 00006 00007 * Redistributions of source code must retain the above copyright notice, this 00008 list of conditions and the following disclaimer. 00009 * Redistributions in binary form must reproduce the above copyright notice, 00010 this list of conditions and the following disclaimer in the documentation 00011 and/or other materials provided with the distribution. 00012 * Neither the name of Intel Corporation nor the names of its contributors may 00013 be used to endorse or promote products derived from this software without 00014 specific prior written permission. 00015 00016 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 00017 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 00018 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 00019 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 00020 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 00021 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 00022 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 00023 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00024 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00025 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00026 00027 ******************************************************************************** 00028 * Content : Eigen bindings to Intel(R) MKL 00029 * Selfadjoint matrix-vector product functionality based on ?SYMV/HEMV. 00030 ******************************************************************************** 00031 */ 00032 00033 #ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H 00034 #define EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H 00035 00036 namespace Eigen { 00037 00038 namespace internal { 00039 00040 /********************************************************************** 00041 * This file implements selfadjoint matrix-vector multiplication using BLAS 00042 **********************************************************************/ 00043 00044 // symv/hemv specialization 00045 00046 template<typename Scalar, typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs> 00047 struct selfadjoint_matrix_vector_product_symv : 00048 selfadjoint_matrix_vector_product<Scalar,Index,StorageOrder,UpLo,ConjugateLhs,ConjugateRhs,BuiltIn> {}; 00049 00050 #define EIGEN_MKL_SYMV_SPECIALIZE(Scalar) \ 00051 template<typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs> \ 00052 struct selfadjoint_matrix_vector_product<Scalar,Index,StorageOrder,UpLo,ConjugateLhs,ConjugateRhs,Specialized> { \ 00053 static void run( \ 00054 Index size, const Scalar* lhs, Index lhsStride, \ 00055 const Scalar* _rhs, Index rhsIncr, Scalar* res, Scalar alpha) { \ 00056 enum {\ 00057 IsColMajor = StorageOrder==ColMajor \ 00058 }; \ 00059 if (IsColMajor == ConjugateLhs) {\ 00060 selfadjoint_matrix_vector_product<Scalar,Index,StorageOrder,UpLo,ConjugateLhs,ConjugateRhs,BuiltIn>::run( \ 00061 size, lhs, lhsStride, _rhs, rhsIncr, res, alpha); \ 00062 } else {\ 00063 selfadjoint_matrix_vector_product_symv<Scalar,Index,StorageOrder,UpLo,ConjugateLhs,ConjugateRhs>::run( \ 00064 size, lhs, lhsStride, _rhs, rhsIncr, res, alpha); \ 00065 }\ 00066 } \ 00067 }; \ 00068 00069 EIGEN_MKL_SYMV_SPECIALIZE(double) 00070 EIGEN_MKL_SYMV_SPECIALIZE(float) 00071 EIGEN_MKL_SYMV_SPECIALIZE(dcomplex) 00072 EIGEN_MKL_SYMV_SPECIALIZE(scomplex) 00073 00074 #define EIGEN_MKL_SYMV_SPECIALIZATION(EIGTYPE,MKLTYPE,MKLFUNC) \ 00075 template<typename Index, int StorageOrder, int UpLo, bool ConjugateLhs, bool ConjugateRhs> \ 00076 struct selfadjoint_matrix_vector_product_symv<EIGTYPE,Index,StorageOrder,UpLo,ConjugateLhs,ConjugateRhs> \ 00077 { \ 00078 typedef Matrix<EIGTYPE,Dynamic,1,ColMajor> SYMVVector;\ 00079 \ 00080 static void run( \ 00081 Index size, const EIGTYPE* lhs, Index lhsStride, \ 00082 const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* res, EIGTYPE alpha) \ 00083 { \ 00084 enum {\ 00085 IsRowMajor = StorageOrder==RowMajor ? 1 : 0, \ 00086 IsLower = UpLo == Lower ? 1 : 0 \ 00087 }; \ 00088 MKL_INT n=size, lda=lhsStride, incx=rhsIncr, incy=1; \ 00089 MKLTYPE alpha_, beta_; \ 00090 const EIGTYPE *x_ptr, myone(1); \ 00091 char uplo=(IsRowMajor) ? (IsLower ? 'U' : 'L') : (IsLower ? 'L' : 'U'); \ 00092 assign_scalar_eig2mkl(alpha_, alpha); \ 00093 assign_scalar_eig2mkl(beta_, myone); \ 00094 SYMVVector x_tmp; \ 00095 if (ConjugateRhs) { \ 00096 Map<const SYMVVector, 0, InnerStride<> > map_x(_rhs,size,1,InnerStride<>(incx)); \ 00097 x_tmp=map_x.conjugate(); \ 00098 x_ptr=x_tmp.data(); \ 00099 incx=1; \ 00100 } else x_ptr=_rhs; \ 00101 MKLFUNC(&uplo, &n, &alpha_, (const MKLTYPE*)lhs, &lda, (const MKLTYPE*)x_ptr, &incx, &beta_, (MKLTYPE*)res, &incy); \ 00102 }\ 00103 }; 00104 00105 EIGEN_MKL_SYMV_SPECIALIZATION(double, double, dsymv) 00106 EIGEN_MKL_SYMV_SPECIALIZATION(float, float, ssymv) 00107 EIGEN_MKL_SYMV_SPECIALIZATION(dcomplex, MKL_Complex16, zhemv) 00108 EIGEN_MKL_SYMV_SPECIALIZATION(scomplex, MKL_Complex8, chemv) 00109 00110 } // end namespace internal 00111 00112 } // end namespace Eigen 00113 00114 #endif // EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H