$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_BLASUTIL_H 00011 #define EIGEN_BLASUTIL_H 00012 00013 // This file contains many lightweight helper classes used to 00014 // implement and control fast level 2 and level 3 BLAS-like routines. 00015 00016 namespace Eigen { 00017 00018 namespace internal { 00019 00020 // forward declarations 00021 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjugateLhs=false, bool ConjugateRhs=false> 00022 struct gebp_kernel; 00023 00024 template<typename Scalar, typename Index, int nr, int StorageOrder, bool Conjugate = false, bool PanelMode=false> 00025 struct gemm_pack_rhs; 00026 00027 template<typename Scalar, typename Index, int Pack1, int Pack2, int StorageOrder, bool Conjugate = false, bool PanelMode = false> 00028 struct gemm_pack_lhs; 00029 00030 template< 00031 typename Index, 00032 typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, 00033 typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs, 00034 int ResStorageOrder> 00035 struct general_matrix_matrix_product; 00036 00037 template<typename Index, typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version=Specialized> 00038 struct general_matrix_vector_product; 00039 00040 00041 template<bool Conjugate> struct conj_if; 00042 00043 template<> struct conj_if<true> { 00044 template<typename T> 00045 inline T operator()(const T& x) { return numext::conj(x); } 00046 template<typename T> 00047 inline T pconj(const T& x) { return internal::pconj(x); } 00048 }; 00049 00050 template<> struct conj_if<false> { 00051 template<typename T> 00052 inline const T& operator()(const T& x) { return x; } 00053 template<typename T> 00054 inline const T& pconj(const T& x) { return x; } 00055 }; 00056 00057 template<typename Scalar> struct conj_helper<Scalar,Scalar,false,false> 00058 { 00059 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); } 00060 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); } 00061 }; 00062 00063 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, false,true> 00064 { 00065 typedef std::complex<RealScalar> Scalar; 00066 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const 00067 { return c + pmul(x,y); } 00068 00069 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const 00070 { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::imag(x)*numext::real(y) - numext::real(x)*numext::imag(y)); } 00071 }; 00072 00073 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false> 00074 { 00075 typedef std::complex<RealScalar> Scalar; 00076 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const 00077 { return c + pmul(x,y); } 00078 00079 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const 00080 { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } 00081 }; 00082 00083 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true> 00084 { 00085 typedef std::complex<RealScalar> Scalar; 00086 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const 00087 { return c + pmul(x,y); } 00088 00089 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const 00090 { return Scalar(numext::real(x)*numext::real(y) - numext::imag(x)*numext::imag(y), - numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } 00091 }; 00092 00093 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false> 00094 { 00095 typedef std::complex<RealScalar> Scalar; 00096 EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const 00097 { return padd(c, pmul(x,y)); } 00098 EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const 00099 { return conj_if<Conj>()(x)*y; } 00100 }; 00101 00102 template<typename RealScalar,bool Conj> struct conj_helper<RealScalar, std::complex<RealScalar>, false,Conj> 00103 { 00104 typedef std::complex<RealScalar> Scalar; 00105 EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const 00106 { return padd(c, pmul(x,y)); } 00107 EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const 00108 { return x*conj_if<Conj>()(y); } 00109 }; 00110 00111 template<typename From,typename To> struct get_factor { 00112 static EIGEN_STRONG_INLINE To run(const From& x) { return x; } 00113 }; 00114 00115 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> { 00116 static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return numext::real(x); } 00117 }; 00118 00119 // Lightweight helper class to access matrix coefficients. 00120 // Yes, this is somehow redundant with Map<>, but this version is much much lighter, 00121 // and so I hope better compilation performance (time and code quality). 00122 template<typename Scalar, typename Index, int StorageOrder> 00123 class blas_data_mapper 00124 { 00125 public: 00126 blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {} 00127 EIGEN_STRONG_INLINE Scalar& operator()(Index i, Index j) 00128 { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; } 00129 protected: 00130 Scalar* EIGEN_RESTRICT m_data; 00131 Index m_stride; 00132 }; 00133 00134 // lightweight helper class to access matrix coefficients (const version) 00135 template<typename Scalar, typename Index, int StorageOrder> 00136 class const_blas_data_mapper 00137 { 00138 public: 00139 const_blas_data_mapper(const Scalar* data, Index stride) : m_data(data), m_stride(stride) {} 00140 EIGEN_STRONG_INLINE const Scalar& operator()(Index i, Index j) const 00141 { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; } 00142 protected: 00143 const Scalar* EIGEN_RESTRICT m_data; 00144 Index m_stride; 00145 }; 00146 00147 00148 /* Helper class to analyze the factors of a Product expression. 00149 * In particular it allows to pop out operator-, scalar multiples, 00150 * and conjugate */ 00151 template<typename XprType> struct blas_traits 00152 { 00153 typedef typename traits<XprType>::Scalar Scalar; 00154 typedef const XprType& ExtractType; 00155 typedef XprType _ExtractType; 00156 enum { 00157 IsComplex = NumTraits<Scalar>::IsComplex, 00158 IsTransposed = false, 00159 NeedToConjugate = false, 00160 HasUsableDirectAccess = ( (int(XprType::Flags)&DirectAccessBit) 00161 && ( bool(XprType::IsVectorAtCompileTime) 00162 || int(inner_stride_at_compile_time<XprType>::ret) == 1) 00163 ) ? 1 : 0 00164 }; 00165 typedef typename conditional<bool(HasUsableDirectAccess), 00166 ExtractType, 00167 typename _ExtractType::PlainObject 00168 >::type DirectLinearAccessType; 00169 static inline ExtractType extract(const XprType& x) { return x; } 00170 static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); } 00171 }; 00172 00173 // pop conjugate 00174 template<typename Scalar, typename NestedXpr> 00175 struct blas_traits<CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> > 00176 : blas_traits<NestedXpr> 00177 { 00178 typedef blas_traits<NestedXpr> Base; 00179 typedef CwiseUnaryOp<scalar_conjugate_op<Scalar>, NestedXpr> XprType; 00180 typedef typename Base::ExtractType ExtractType; 00181 00182 enum { 00183 IsComplex = NumTraits<Scalar>::IsComplex, 00184 NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex 00185 }; 00186 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } 00187 static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); } 00188 }; 00189 00190 // pop scalar multiple 00191 template<typename Scalar, typename NestedXpr> 00192 struct blas_traits<CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> > 00193 : blas_traits<NestedXpr> 00194 { 00195 typedef blas_traits<NestedXpr> Base; 00196 typedef CwiseUnaryOp<scalar_multiple_op<Scalar>, NestedXpr> XprType; 00197 typedef typename Base::ExtractType ExtractType; 00198 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } 00199 static inline Scalar extractScalarFactor(const XprType& x) 00200 { return x.functor().m_other * Base::extractScalarFactor(x.nestedExpression()); } 00201 }; 00202 00203 // pop opposite 00204 template<typename Scalar, typename NestedXpr> 00205 struct blas_traits<CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> > 00206 : blas_traits<NestedXpr> 00207 { 00208 typedef blas_traits<NestedXpr> Base; 00209 typedef CwiseUnaryOp<scalar_opposite_op<Scalar>, NestedXpr> XprType; 00210 typedef typename Base::ExtractType ExtractType; 00211 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } 00212 static inline Scalar extractScalarFactor(const XprType& x) 00213 { return - Base::extractScalarFactor(x.nestedExpression()); } 00214 }; 00215 00216 // pop/push transpose 00217 template<typename NestedXpr> 00218 struct blas_traits<Transpose<NestedXpr> > 00219 : blas_traits<NestedXpr> 00220 { 00221 typedef typename NestedXpr::Scalar Scalar; 00222 typedef blas_traits<NestedXpr> Base; 00223 typedef Transpose<NestedXpr> XprType; 00224 typedef Transpose<const typename Base::_ExtractType> ExtractType; // const to get rid of a compile error; anyway blas traits are only used on the RHS 00225 typedef Transpose<const typename Base::_ExtractType> _ExtractType; 00226 typedef typename conditional<bool(Base::HasUsableDirectAccess), 00227 ExtractType, 00228 typename ExtractType::PlainObject 00229 >::type DirectLinearAccessType; 00230 enum { 00231 IsTransposed = Base::IsTransposed ? 0 : 1 00232 }; 00233 static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } 00234 static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); } 00235 }; 00236 00237 template<typename T> 00238 struct blas_traits<const T> 00239 : blas_traits<T> 00240 {}; 00241 00242 template<typename T, bool HasUsableDirectAccess=blas_traits<T>::HasUsableDirectAccess> 00243 struct extract_data_selector { 00244 static const typename T::Scalar* run(const T& m) 00245 { 00246 return blas_traits<T>::extract(m).data(); 00247 } 00248 }; 00249 00250 template<typename T> 00251 struct extract_data_selector<T,false> { 00252 static typename T::Scalar* run(const T&) { return 0; } 00253 }; 00254 00255 template<typename T> const typename T::Scalar* extract_data(const T& m) 00256 { 00257 return extract_data_selector<T>::run(m); 00258 } 00259 00260 } // end namespace internal 00261 00262 } // end namespace Eigen 00263 00264 #endif // EIGEN_BLASUTIL_H