Main MRPT website > C++ reference
MRPT logo

BlasUtil.h

Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_BLASUTIL_H
00026 #define EIGEN_BLASUTIL_H
00027 
00028 // This file contains many lightweight helper classes used to
00029 // implement and control fast level 2 and level 3 BLAS-like routines.
00030 
00031 namespace internal {
00032 
00033 // forward declarations
00034 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjugateLhs=false, bool ConjugateRhs=false>
00035 struct gebp_kernel;
00036 
00037 template<typename Scalar, typename Index, int nr, int StorageOrder, bool Conjugate = false, bool PanelMode=false>
00038 struct gemm_pack_rhs;
00039 
00040 template<typename Scalar, typename Index, int Pack1, int Pack2, int StorageOrder, bool Conjugate = false, bool PanelMode = false>
00041 struct gemm_pack_lhs;
00042 
00043 template<
00044   typename Index,
00045   typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00046   typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs,
00047   int ResStorageOrder>
00048 struct general_matrix_matrix_product;
00049 
00050 template<typename Index, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs>
00051 struct general_matrix_vector_product;
00052 
00053 
00054 template<bool Conjugate> struct conj_if;
00055 
00056 template<> struct conj_if<true> {
00057   template<typename T>
00058   inline T operator()(const T& x) { return conj(x); }
00059 };
00060 
00061 template<> struct conj_if<false> {
00062   template<typename T>
00063   inline const T& operator()(const T& x) { return x; }
00064 };
00065 
00066 template<typename Scalar> struct conj_helper<Scalar,Scalar,false,false>
00067 {
00068   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); }
00069   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); }
00070 };
00071 
00072 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, false,true>
00073 {
00074   typedef std::complex<RealScalar> Scalar;
00075   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
00076   { return c + pmul(x,y); }
00077 
00078   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
00079   { return Scalar(real(x)*real(y) + imag(x)*imag(y), imag(x)*real(y) - real(x)*imag(y)); }
00080 };
00081 
00082 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false>
00083 {
00084   typedef std::complex<RealScalar> Scalar;
00085   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
00086   { return c + pmul(x,y); }
00087 
00088   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
00089   { return Scalar(real(x)*real(y) + imag(x)*imag(y), real(x)*imag(y) - imag(x)*real(y)); }
00090 };
00091 
00092 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true>
00093 {
00094   typedef std::complex<RealScalar> Scalar;
00095   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const
00096   { return c + pmul(x,y); }
00097 
00098   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
00099   { return Scalar(real(x)*real(y) - imag(x)*imag(y), - real(x)*imag(y) - imag(x)*real(y)); }
00100 };
00101 
00102 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false>
00103 {
00104   typedef std::complex<RealScalar> Scalar;
00105   EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const
00106   { return padd(c, pmul(x,y)); }
00107   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const
00108   { return conj_if<Conj>()(x)*y; }
00109 };
00110 
00111 template<typename RealScalar,bool Conj> struct conj_helper<RealScalar, std::complex<RealScalar>, false,Conj>
00112 {
00113   typedef std::complex<RealScalar> Scalar;
00114   EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const
00115   { return padd(c, pmul(x,y)); }
00116   EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const
00117   { return x*conj_if<Conj>()(y); }
00118 };
00119 
00120 template<typename From,typename To> struct get_factor {
00121   EIGEN_STRONG_INLINE static To run(const From& x) { return x; }
00122 };
00123 
00124 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> {
00125   EIGEN_STRONG_INLINE static typename NumTraits<Scalar>::Real run(const Scalar& x) { return real(x); }
00126 };
00127 
00128 // Lightweight helper class to access matrix coefficients.
00129 // Yes, this is somehow redundant with Map<>, but this version is much much lighter,
00130 // and so I hope better compilation performance (time and code quality).
00131 template<typename Scalar, typename Index, int StorageOrder>
00132 class blas_data_mapper
00133 {
00134   public:
00135     blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
00136     EIGEN_STRONG_INLINE Scalar& operator()(Index i, Index j)
00137     { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
00138   protected:
00139     Scalar* EIGEN_RESTRICT m_data;
00140     Index m_stride;
00141 };
00142 
00143 // lightweight helper class to access matrix coefficients (const version)
00144 template<typename Scalar, typename Index, int StorageOrder>
00145 class const_blas_data_mapper
00146 {
00147   public:
00148     const_blas_data_mapper(const Scalar* data, Index stride) : m_data(data), m_stride(stride) {}
00149     EIGEN_STRONG_INLINE const Scalar& operator()(Index i, Index j) const
00150     { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; }
00151   protected:
00152     const Scalar* EIGEN_RESTRICT m_data;
00153     Index m_stride;
00154 };
00155 
00156 
00157 /* Helper class to analyze the factors of a Product expression.
00158  * In particular it allows to pop out operator-, scalar multiples,
00159  * and conjugate */
00160 template<typename XprType> struct blas_traits
00161 {
00162   typedef typename traits<XprType>::Scalar Scalar;
00163   typedef const XprType& ExtractType;
00164   typedef XprType _ExtractType;
00165   enum {
00166     IsComplex = NumTraits<Scalar>::IsComplex,
00167     IsTransposed = false,
00168     NeedToConjugate = false,
00169     HasUsableDirectAccess = (    (int(XprType::Flags)&DirectAccessBit)
00170                      && (  /* Uncomment this when the low-level matrix-vector product functions support strided vectors
00171                            bool(XprType::IsVectorAtCompileTime)
00172                          ||  */
00173                            int(inner_stride_at_compile_time<XprType>::ret) == 1)
00174                    ) ?  1 : 0
00175   };
00176   typedef typename conditional<bool(HasUsableDirectAccess),
00177     ExtractType,
00178     typename _ExtractType::PlainObject
00179     >::type DirectLinearAccessType;
00180   static inline const ExtractType extract(const XprType& x) { return x; }
00181   static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
00182 };
00183 
00184 // pop conjugate
00185 template<typename Scalar, typename NestedXpr>
00186 struct blas_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> >
00187  : blas_traits<NestedXpr>
00188 {
00189   typedef blas_traits<NestedXpr> Base;
00190   typedef CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> XprType;
00191   typedef typename Base::ExtractType ExtractType;
00192 
00193   enum {
00194     IsComplex = NumTraits<Scalar>::IsComplex,
00195     NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex
00196   };
00197   static inline const ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
00198   static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); }
00199 };
00200 
00201 // pop scalar multiple
00202 template<typename Scalar, typename NestedXpr>
00203 struct blas_traits<CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> >
00204  : blas_traits<NestedXpr>
00205 {
00206   typedef blas_traits<NestedXpr> Base;
00207   typedef CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> XprType;
00208   typedef typename Base::ExtractType ExtractType;
00209   static inline const ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
00210   static inline Scalar extractScalarFactor(const XprType& x)
00211   { return x.functor().m_other * Base::extractScalarFactor(x.nestedExpression()); }
00212 };
00213 
00214 // pop opposite
00215 template<typename Scalar, typename NestedXpr>
00216 struct blas_traits<CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> >
00217  : blas_traits<NestedXpr>
00218 {
00219   typedef blas_traits<NestedXpr> Base;
00220   typedef CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> XprType;
00221   typedef typename Base::ExtractType ExtractType;
00222   static inline const ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
00223   static inline Scalar extractScalarFactor(const XprType& x)
00224   { return - Base::extractScalarFactor(x.nestedExpression()); }
00225 };
00226 
00227 // pop/push transpose
00228 template<typename NestedXpr>
00229 struct blas_traits<Transpose<NestedXpr> >
00230  : blas_traits<NestedXpr>
00231 {
00232   typedef typename NestedXpr::Scalar Scalar;
00233   typedef blas_traits<NestedXpr> Base;
00234   typedef Transpose<NestedXpr> XprType;
00235   typedef Transpose<typename Base::_ExtractType>  ExtractType;
00236   typedef Transpose<typename Base::_ExtractType> _ExtractType;
00237   typedef typename conditional<bool(Base::HasUsableDirectAccess),
00238     ExtractType,
00239     typename ExtractType::PlainObject
00240     >::type DirectLinearAccessType;
00241   enum {
00242     IsTransposed = Base::IsTransposed ? 0 : 1
00243   };
00244   static inline const ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); }
00245   static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); }
00246 };
00247 
00248 template<typename T>
00249 struct blas_traits<const T>
00250      : blas_traits<T>
00251 {};
00252 
00253 template<typename T, bool HasUsableDirectAccess=blas_traits<T>::HasUsableDirectAccess>
00254 struct extract_data_selector {
00255   static const typename T::Scalar* run(const T& m)
00256   {
00257     return const_cast<typename T::Scalar*>(&blas_traits<T>::extract(m).coeffRef(0,0)); // FIXME this should be .data()
00258   }
00259 };
00260 
00261 template<typename T>
00262 struct extract_data_selector<T,false> {
00263   static typename T::Scalar* run(const T&) { return 0; }
00264 };
00265 
00266 template<typename T> const typename T::Scalar* extract_data(const T& m)
00267 {
00268   return extract_data_selector<T>::run(m);
00269 }
00270 
00271 } // end namespace internal
00272 
00273 #endif // EIGEN_BLASUTIL_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011