Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
TriangularMatrixVector_MKL.h
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 * Triangular matrix-vector product functionality based on ?TRMV. 00030 ******************************************************************************** 00031 */ 00032 00033 #ifndef EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H 00034 #define EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H 00035 00036 namespace Eigen { 00037 00038 namespace internal { 00039 00040 /********************************************************************** 00041 * This file implements triangular matrix-vector multiplication using BLAS 00042 **********************************************************************/ 00043 00044 // trmv/hemv specialization 00045 00046 template<typename Index, int Mode, typename LhsScalar, bool ConjLhs, typename RhsScalar, bool ConjRhs, int StorageOrder> 00047 struct triangular_matrix_vector_product_trmv : 00048 triangular_matrix_vector_product<Index,Mode,LhsScalar,ConjLhs,RhsScalar,ConjRhs,StorageOrder,BuiltIn> {}; 00049 00050 #define EIGEN_MKL_TRMV_SPECIALIZE(Scalar) \ 00051 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \ 00052 struct triangular_matrix_vector_product<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,ColMajor,Specialized> { \ 00053 static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ 00054 const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ 00055 triangular_matrix_vector_product_trmv<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,ColMajor>::run( \ 00056 _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ 00057 } \ 00058 }; \ 00059 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \ 00060 struct triangular_matrix_vector_product<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,RowMajor,Specialized> { \ 00061 static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ 00062 const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ 00063 triangular_matrix_vector_product_trmv<Index,Mode,Scalar,ConjLhs,Scalar,ConjRhs,RowMajor>::run( \ 00064 _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ 00065 } \ 00066 }; 00067 00068 EIGEN_MKL_TRMV_SPECIALIZE(double) 00069 EIGEN_MKL_TRMV_SPECIALIZE(float) 00070 EIGEN_MKL_TRMV_SPECIALIZE(dcomplex) 00071 EIGEN_MKL_TRMV_SPECIALIZE(scomplex) 00072 00073 // implements col-major: res += alpha * op(triangular) * vector 00074 #define EIGEN_MKL_TRMV_CM(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 00075 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \ 00076 struct triangular_matrix_vector_product_trmv<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,ColMajor> { \ 00077 enum { \ 00078 IsLower = (Mode&Lower) == Lower, \ 00079 SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ 00080 IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ 00081 IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ 00082 LowUp = IsLower ? Lower : Upper \ 00083 }; \ 00084 static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ 00085 const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ 00086 { \ 00087 if (ConjLhs || IsZeroDiag) { \ 00088 triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,ColMajor,BuiltIn>::run( \ 00089 _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ 00090 return; \ 00091 }\ 00092 Index size = (std::min)(_rows,_cols); \ 00093 Index rows = IsLower ? _rows : size; \ 00094 Index cols = IsLower ? size : _cols; \ 00095 \ 00096 typedef VectorX##EIGPREFIX VectorRhs; \ 00097 EIGTYPE *x, *y;\ 00098 \ 00099 /* Set x*/ \ 00100 Map<const VectorRhs, 0, InnerStride<> > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ 00101 VectorRhs x_tmp; \ 00102 if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ 00103 x = x_tmp.data(); \ 00104 \ 00105 /* Square part handling */\ 00106 \ 00107 char trans, uplo, diag; \ 00108 MKL_INT m, n, lda, incx, incy; \ 00109 EIGTYPE const *a; \ 00110 MKLTYPE alpha_, beta_; \ 00111 assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(alpha_, alpha); \ 00112 assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(beta_, EIGTYPE(1)); \ 00113 \ 00114 /* Set m, n */ \ 00115 n = (MKL_INT)size; \ 00116 lda = lhsStride; \ 00117 incx = 1; \ 00118 incy = resIncr; \ 00119 \ 00120 /* Set uplo, trans and diag*/ \ 00121 trans = 'N'; \ 00122 uplo = IsLower ? 'L' : 'U'; \ 00123 diag = IsUnitDiag ? 'U' : 'N'; \ 00124 \ 00125 /* call ?TRMV*/ \ 00126 MKLPREFIX##trmv(&uplo, &trans, &diag, &n, (const MKLTYPE*)_lhs, &lda, (MKLTYPE*)x, &incx); \ 00127 \ 00128 /* Add op(a_tr)rhs into res*/ \ 00129 MKLPREFIX##axpy(&n, &alpha_,(const MKLTYPE*)x, &incx, (MKLTYPE*)_res, &incy); \ 00130 /* Non-square case - doesn't fit to MKL ?TRMV. Fall to default triangular product*/ \ 00131 if (size<(std::max)(rows,cols)) { \ 00132 typedef Matrix<EIGTYPE, Dynamic, Dynamic> MatrixLhs; \ 00133 if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ 00134 x = x_tmp.data(); \ 00135 if (size<rows) { \ 00136 y = _res + size*resIncr; \ 00137 a = _lhs + size; \ 00138 m = rows-size; \ 00139 n = size; \ 00140 } \ 00141 else { \ 00142 x += size; \ 00143 y = _res; \ 00144 a = _lhs + size*lda; \ 00145 m = size; \ 00146 n = cols-size; \ 00147 } \ 00148 MKLPREFIX##gemv(&trans, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)x, &incx, &beta_, (MKLTYPE*)y, &incy); \ 00149 } \ 00150 } \ 00151 }; 00152 00153 EIGEN_MKL_TRMV_CM(double, double, d, d) 00154 EIGEN_MKL_TRMV_CM(dcomplex, MKL_Complex16, cd, z) 00155 EIGEN_MKL_TRMV_CM(float, float, f, s) 00156 EIGEN_MKL_TRMV_CM(scomplex, MKL_Complex8, cf, c) 00157 00158 // implements row-major: res += alpha * op(triangular) * vector 00159 #define EIGEN_MKL_TRMV_RM(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 00160 template<typename Index, int Mode, bool ConjLhs, bool ConjRhs> \ 00161 struct triangular_matrix_vector_product_trmv<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,RowMajor> { \ 00162 enum { \ 00163 IsLower = (Mode&Lower) == Lower, \ 00164 SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ 00165 IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ 00166 IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ 00167 LowUp = IsLower ? Lower : Upper \ 00168 }; \ 00169 static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ 00170 const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ 00171 { \ 00172 if (IsZeroDiag) { \ 00173 triangular_matrix_vector_product<Index,Mode,EIGTYPE,ConjLhs,EIGTYPE,ConjRhs,RowMajor,BuiltIn>::run( \ 00174 _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ 00175 return; \ 00176 }\ 00177 Index size = (std::min)(_rows,_cols); \ 00178 Index rows = IsLower ? _rows : size; \ 00179 Index cols = IsLower ? size : _cols; \ 00180 \ 00181 typedef VectorX##EIGPREFIX VectorRhs; \ 00182 EIGTYPE *x, *y;\ 00183 \ 00184 /* Set x*/ \ 00185 Map<const VectorRhs, 0, InnerStride<> > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ 00186 VectorRhs x_tmp; \ 00187 if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ 00188 x = x_tmp.data(); \ 00189 \ 00190 /* Square part handling */\ 00191 \ 00192 char trans, uplo, diag; \ 00193 MKL_INT m, n, lda, incx, incy; \ 00194 EIGTYPE const *a; \ 00195 MKLTYPE alpha_, beta_; \ 00196 assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(alpha_, alpha); \ 00197 assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(beta_, EIGTYPE(1)); \ 00198 \ 00199 /* Set m, n */ \ 00200 n = (MKL_INT)size; \ 00201 lda = lhsStride; \ 00202 incx = 1; \ 00203 incy = resIncr; \ 00204 \ 00205 /* Set uplo, trans and diag*/ \ 00206 trans = ConjLhs ? 'C' : 'T'; \ 00207 uplo = IsLower ? 'U' : 'L'; \ 00208 diag = IsUnitDiag ? 'U' : 'N'; \ 00209 \ 00210 /* call ?TRMV*/ \ 00211 MKLPREFIX##trmv(&uplo, &trans, &diag, &n, (const MKLTYPE*)_lhs, &lda, (MKLTYPE*)x, &incx); \ 00212 \ 00213 /* Add op(a_tr)rhs into res*/ \ 00214 MKLPREFIX##axpy(&n, &alpha_,(const MKLTYPE*)x, &incx, (MKLTYPE*)_res, &incy); \ 00215 /* Non-square case - doesn't fit to MKL ?TRMV. Fall to default triangular product*/ \ 00216 if (size<(std::max)(rows,cols)) { \ 00217 typedef Matrix<EIGTYPE, Dynamic, Dynamic> MatrixLhs; \ 00218 if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ 00219 x = x_tmp.data(); \ 00220 if (size<rows) { \ 00221 y = _res + size*resIncr; \ 00222 a = _lhs + size*lda; \ 00223 m = rows-size; \ 00224 n = size; \ 00225 } \ 00226 else { \ 00227 x += size; \ 00228 y = _res; \ 00229 a = _lhs + size; \ 00230 m = size; \ 00231 n = cols-size; \ 00232 } \ 00233 MKLPREFIX##gemv(&trans, &n, &m, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)x, &incx, &beta_, (MKLTYPE*)y, &incy); \ 00234 } \ 00235 } \ 00236 }; 00237 00238 EIGEN_MKL_TRMV_RM(double, double, d, d) 00239 EIGEN_MKL_TRMV_RM(dcomplex, MKL_Complex16, cd, z) 00240 EIGEN_MKL_TRMV_RM(float, float, f, s) 00241 EIGEN_MKL_TRMV_RM(scomplex, MKL_Complex8, cf, c) 00242 00243 } // end namespase internal 00244 00245 } // end namespace Eigen 00246 00247 #endif // EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H
Generated on Tue Jul 12 2022 17:47:01 by 1.7.2