turtlebot_exploration_3d: SelfadjointMatrixVector

Go to the documentation of this file.
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