diff --git a/Eigen/Core b/Eigen/Core index ea871dca3..d6a72765e 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -176,6 +176,7 @@ namespace Eigen { #include "src/Core/IO.h" #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" +#include "src/Core/util/BlasUtil.h" #include "src/Core/Product.h" #include "src/Core/products/GeneralMatrixMatrix.h" #include "src/Core/products/GeneralMatrixVector.h" diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 93e318035..78cb88c33 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -577,22 +577,6 @@ struct ei_product_packet_impl -static void ei_cache_friendly_product( - int _rows, int _cols, int depth, - bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride, - bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride, - bool resRowMajor, Scalar* res, int resStride, - Scalar alpha); - -template -static void ei_cache_friendly_product_colmajor_times_vector( - int size, const Scalar* lhs, int lhsStride, const RhsType& rhs, Scalar* res, Scalar alpha); - -template -static void ei_cache_friendly_product_rowmajor_times_vector( - const Scalar* lhs, int lhsStride, const Scalar* rhs, int rhsSize, ResType& res, Scalar alpha); - // This helper class aims to determine which optimized product to call, // and how to call it. We have to distinghish three major cases: // 1 - matrix-matrix @@ -926,16 +910,18 @@ inline void Product::_cacheFriendlyEvalAndAdd(DestDerived& typedef typename ei_unref::type _RhsCopy; LhsCopy lhs(actualLhs); RhsCopy rhs(actualRhs); - ei_cache_friendly_product - ( - rows(), cols(), lhs.cols(), - _LhsCopy::Flags&RowMajorBit, (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(), - _RhsCopy::Flags&RowMajorBit, (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(), - Flags&RowMajorBit, (Scalar*)&(res.coeffRef(0,0)), res.stride(), - actualAlpha - ); + + ei_general_matrix_matrix_product< + Scalar, + (_LhsCopy::Flags&RowMajorBit)?RowMajor:ColMajor, bool(LhsProductTraits::NeedToConjugate), + (_RhsCopy::Flags&RowMajorBit)?RowMajor:ColMajor, bool(RhsProductTraits::NeedToConjugate), + (DestDerived::Flags&RowMajorBit)?RowMajor:ColMajor> + ::run( + rows(), cols(), lhs.cols(), + (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(), + (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(), + (Scalar*)&(res.coeffRef(0,0)), res.stride(), + actualAlpha); } #endif // EIGEN_PRODUCT_H diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 30fa4aa66..68949499a 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -25,145 +25,63 @@ #ifndef EIGEN_GENERAL_MATRIX_MATRIX_H #define EIGEN_GENERAL_MATRIX_MATRIX_H -template -struct ei_gebp_kernel; - -template -struct ei_gemm_pack_rhs; - -template -struct ei_gemm_pack_lhs; - -template -struct ei_L2_block_traits { - enum {width = 8 * ei_meta_sqrt::ret }; -}; - -template struct ei_conj_helper; - -template<> struct ei_conj_helper -{ - template - EIGEN_STRONG_INLINE T pmadd(const T& x, const T& y, const T& c) const { return ei_pmadd(x,y,c); } - template - EIGEN_STRONG_INLINE T pmul(const T& x, const T& y) const { return ei_pmul(x,y); } -}; - -template<> struct ei_conj_helper -{ - template std::complex - pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const - { return c + pmul(x,y); } - - template std::complex pmul(const std::complex& x, const std::complex& y) const - { return std::complex(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_imag(x)*ei_real(y) - ei_real(x)*ei_imag(y)); } -}; - -template<> struct ei_conj_helper -{ - template std::complex - pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const - { return c + pmul(x,y); } - - template std::complex pmul(const std::complex& x, const std::complex& y) const - { return std::complex(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); } -}; - -template<> struct ei_conj_helper -{ - template std::complex - pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const - { return c + pmul(x,y); } - - template std::complex pmul(const std::complex& x, const std::complex& y) const - { return std::complex(ei_real(x)*ei_real(y) - ei_imag(x)*ei_imag(y), - ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); } -}; - #ifndef EIGEN_EXTERN_INSTANTIATIONS -/** \warning you should never call this function directly, - * this is because the ConjugateLhs/ConjugateRhs have to - * be flipped is resRowMajor==true */ -template -static void ei_cache_friendly_product( - int _rows, int _cols, int depth, - bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride, - bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride, - bool resRowMajor, Scalar* res, int resStride, +template< + typename Scalar, + int LhsStorageOrder, bool ConjugateLhs, + int RhsStorageOrder, bool ConjugateRhs> +struct ei_general_matrix_matrix_product +{ + static EIGEN_STRONG_INLINE void run( + int rows, int cols, int depth, + const Scalar* lhs, int lhsStride, + const Scalar* rhs, int rhsStride, + Scalar* res, int resStride, + Scalar alpha) + { + // transpose the product such that the result is column major + ei_general_matrix_matrix_product + ::run(cols,rows,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha); + } +}; + +template< + typename Scalar, + int LhsStorageOrder, bool ConjugateLhs, + int RhsStorageOrder, bool ConjugateRhs> +struct ei_general_matrix_matrix_product +{ +static void run(int rows, int cols, int depth, + const Scalar* _lhs, int lhsStride, + const Scalar* _rhs, int rhsStride, + Scalar* res, int resStride, Scalar alpha) { - const Scalar* EIGEN_RESTRICT lhs; - const Scalar* EIGEN_RESTRICT rhs; - int lhsStride, rhsStride, rows, cols; - bool lhsRowMajor; + ei_const_blas_data_mapper lhs(_lhs,lhsStride); + ei_const_blas_data_mapper rhs(_rhs,rhsStride); - ei_conj_helper cj; if (ConjugateRhs) alpha = ei_conj(alpha); - bool hasAlpha = alpha != Scalar(1); - - if (resRowMajor) - { -// return ei_cache_friendly_product(_cols,_rows,depth, -// !_rhsRowMajor, _rhs, _rhsStride, -// !_lhsRowMajor, _lhs, _lhsStride, -// false, res, resStride, -// alpha); - - lhs = _rhs; - rhs = _lhs; - lhsStride = _rhsStride; - rhsStride = _lhsStride; - cols = _rows; - rows = _cols; - lhsRowMajor = !_rhsRowMajor; - ei_assert(_lhsRowMajor); - } - else - { - lhs = _lhs; - rhs = _rhs; - lhsStride = _lhsStride; - rhsStride = _rhsStride; - rows = _rows; - cols = _cols; - lhsRowMajor = _lhsRowMajor; - ei_assert(!_rhsRowMajor); - } typedef typename ei_packet_traits::type PacketType; + typedef ei_product_blocking_traits Blocking; - enum { - PacketSize = sizeof(PacketType)/sizeof(Scalar), - #if (defined __i386__) - HalfRegisterCount = 4, - #else - HalfRegisterCount = 8, - #endif - - // register block size along the N direction - nr = HalfRegisterCount/2, - - // register block size along the M direction - mr = 2 * PacketSize, - - // max cache block size along the K direction - Max_kc = ei_L2_block_traits::width, - - // max cache block size along the M direction - Max_mc = 2*Max_kc - }; - - int kc = std::min(Max_kc,depth); // cache block size along the K direction - int mc = std::min(Max_mc,rows); // cache block size along the M direction + int kc = std::min(Blocking::Max_kc,depth); // cache block size along the K direction + int mc = std::min(Blocking::Max_mc,rows); // cache block size along the M direction Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); - Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*PacketSize); + Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*Blocking::PacketSize); // number of columns which can be processed by packet of nr columns - int packet_cols = (cols/nr)*nr; + int packet_cols = (cols/Blocking::nr) * Blocking::nr; - // GEMM_VAR1 + // => GEMM_VAR1 for(int k2=0; k2()(blockB, rhs, rhsStride, hasAlpha, alpha, actual_kc, packet_cols, k2, cols); + ei_gemm_pack_rhs()(blockB, &rhs(k2,0), rhsStride, alpha, actual_kc, packet_cols, cols); // => GEPP_VAR1 for(int i2=0; i2()(blockA, &lhs(i2,k2), lhsStride, actual_kc, actual_mc); - ei_gemm_pack_lhs()(blockA, lhs, lhsStride, lhsRowMajor, actual_kc, actual_mc, k2, i2); - - ei_gebp_kernel >() + ei_gebp_kernel >() (res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); } } ei_aligned_stack_delete(Scalar, blockA, kc*mc); - ei_aligned_stack_delete(Scalar, blockB, kc*cols*PacketSize); + ei_aligned_stack_delete(Scalar, blockB, kc*cols*Blocking::PacketSize); } +}; + // optimized GEneral packed Block * packed Panel product kernel template struct ei_gebp_kernel @@ -356,8 +276,7 @@ struct ei_gebp_kernel if(nr==4) res[(j2+3)*resStride + i2 + i] += C3; } } - // remaining rhs/res columns ( do the same but with nr==1 for(int j2=packet_cols; j2 +template struct ei_gemm_pack_lhs { - void operator()(Scalar* blockA, const Scalar* lhs, int lhsStride, bool lhsRowMajor, int actual_kc, int actual_mc, int k2, int i2) + void operator()(Scalar* blockA, const EIGEN_RESTRICT Scalar* _lhs, int lhsStride, int actual_kc, int actual_mc) { + ei_const_blas_data_mapper lhs(_lhs,lhsStride); int count = 0; const int peeled_mc = (actual_mc/mr)*mr; - if (lhsRowMajor) + for(int i=0; i struct ei_gemm_pack_rhs { - void operator()(Scalar* blockB, const Scalar* rhs, int rhsStride, bool hasAlpha, Scalar alpha, int actual_kc, int packet_cols, int k2, int cols) + void operator()(Scalar* blockB, const Scalar* rhs, int rhsStride, Scalar alpha, int actual_kc, int packet_cols, int cols) { + bool hasAlpha = alpha != Scalar(1); int count = 0; for(int j2=0; j2 +// pack a selfadjoint block diagonal for use with the gebp_kernel +template struct ei_symm_pack_lhs { - void operator()(Scalar* blockA, const Scalar* lhs, int lhsStride, bool lhsRowMajor, int actual_kc, int actual_mc, int k2, int i2) + void operator()(Scalar* blockA, const Scalar* _lhs, int lhsStride, int actual_kc, int actual_mc) { + ei_const_blas_data_mapper lhs(_lhs,lhsStride); int count = 0; const int peeled_mc = (actual_mc/mr)*mr; - if (lhsRowMajor) + for(int i=0; i static EIGEN_DONT_INLINE void ei_product_selfadjoint_matrix( int size, - const Scalar* lhs, int lhsStride, - const Scalar* rhs, int rhsStride, bool rhsRowMajor, int cols, + const Scalar* _lhs, int lhsStride, + const Scalar* _rhs, int rhsStride, bool rhsRowMajor, int cols, Scalar* res, int resStride, Scalar alpha) { typedef typename ei_packet_traits::type Packet; - ei_conj_helper cj; + ei_const_blas_data_mapper lhs(_lhs,lhsStride); + ei_const_blas_data_mapper rhs(_rhs,rhsStride); + if (ConjugateRhs) alpha = ei_conj(alpha); - bool hasAlpha = alpha != Scalar(1); typedef typename ei_packet_traits::type PacketType; const bool lhsRowMajor = (StorageOrder==RowMajor); - enum { - PacketSize = sizeof(PacketType)/sizeof(Scalar), - #if (defined __i386__) - HalfRegisterCount = 4, - #else - HalfRegisterCount = 8, - #endif + typedef ei_product_blocking_traits Blocking; - // register block size along the N direction - nr = HalfRegisterCount/2, - - // register block size along the M direction - mr = 2 * PacketSize, - - // max cache block size along the K direction - Max_kc = ei_L2_block_traits::width, - - // max cache block size along the M direction - Max_mc = 2*Max_kc - }; - - int kc = std::min(Max_kc,size); // cache block size along the K direction - int mc = std::min(Max_mc,size); // cache block size along the M direction + int kc = std::min(Blocking::Max_kc,size); // cache block size along the K direction + int mc = std::min(Blocking::Max_mc,size); // cache block size along the M direction Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); - Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*PacketSize); + Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*Blocking::PacketSize); // number of columns which can be processed by packet of nr columns - int packet_cols = (cols/nr)*nr; + int packet_cols = (cols/Blocking::nr)*Blocking::nr; + + ei_gebp_kernel > gebp_kernel; for(int k2=0; k2()(blockB, rhs, rhsStride, hasAlpha, alpha, actual_kc, packet_cols, k2, cols); + ei_gemm_pack_rhs() + (blockB, &rhs(k2,0), rhsStride, alpha, actual_kc, packet_cols, cols); // the select lhs's panel has to be split in three different parts: // 1 - the transposed panel above the diagonal block => transposed packed copy @@ -153,32 +122,31 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_matrix( { const int actual_mc = std::min(i2+mc,k2)-i2; // transposed packed copy - ei_gemm_pack_lhs()(blockA, lhs, lhsStride, !lhsRowMajor, actual_kc, actual_mc, k2, i2); + ei_gemm_pack_lhs() + (blockA, &lhs(k2,i2), lhsStride, actual_kc, actual_mc); - ei_gebp_kernel >() - (res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); + gebp_kernel(res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); } // the block diagonal { const int actual_mc = std::min(k2+kc,size)-k2; // symmetric packed copy - ei_symm_pack_lhs()(blockA, lhs, lhsStride, lhsRowMajor, actual_kc, actual_mc, k2, k2); - ei_gebp_kernel >() - (res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, k2, cols); + ei_symm_pack_lhs() + (blockA, &lhs(k2,k2), lhsStride, actual_kc, actual_mc); + gebp_kernel(res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, k2, cols); } for(int i2=k2+kc; i2()(blockA, lhs, lhsStride, lhsRowMajor, actual_kc, actual_mc, k2, i2); - ei_gebp_kernel >() - (res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); + ei_gemm_pack_lhs() + (blockA, &lhs(i2,k2), lhsStride, actual_kc, actual_mc); + gebp_kernel(res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); } } ei_aligned_stack_delete(Scalar, blockA, kc*mc); - ei_aligned_stack_delete(Scalar, blockB, kc*cols*PacketSize); + ei_aligned_stack_delete(Scalar, blockB, kc*cols*Blocking::PacketSize); } - #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_H diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h new file mode 100644 index 000000000..6e4b21e6a --- /dev/null +++ b/Eigen/src/Core/util/BlasUtil.h @@ -0,0 +1,157 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_BLASUTIL_H +#define EIGEN_BLASUTIL_H + +// This file contains many lightweight helper classes used to +// implement and control fast level 2 and level 3 BLAS-like routines. + +// forward declarations +template +struct ei_gebp_kernel; + +template +struct ei_gemm_pack_rhs; + +template +struct ei_gemm_pack_lhs; + +template< + typename Scalar, + int LhsStorageOrder, bool ConjugateLhs, + int RhsStorageOrder, bool ConjugateRhs, + int ResStorageOrder> +struct ei_general_matrix_matrix_product; + +template +static void ei_cache_friendly_product_colmajor_times_vector( + int size, const Scalar* lhs, int lhsStride, const RhsType& rhs, Scalar* res, Scalar alpha); + +template +static void ei_cache_friendly_product_rowmajor_times_vector( + const Scalar* lhs, int lhsStride, const Scalar* rhs, int rhsSize, ResType& res, Scalar alpha); + +// Provides scalar/packet-wise product and product with accumulation +// with optional conjugation of the arguments. +template struct ei_conj_helper; + +template<> struct ei_conj_helper +{ + template + EIGEN_STRONG_INLINE T pmadd(const T& x, const T& y, const T& c) const { return ei_pmadd(x,y,c); } + template + EIGEN_STRONG_INLINE T pmul(const T& x, const T& y) const { return ei_pmul(x,y); } +}; + +template<> struct ei_conj_helper +{ + template std::complex + pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const + { return c + pmul(x,y); } + + template std::complex pmul(const std::complex& x, const std::complex& y) const + { return std::complex(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_imag(x)*ei_real(y) - ei_real(x)*ei_imag(y)); } +}; + +template<> struct ei_conj_helper +{ + template std::complex + pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const + { return c + pmul(x,y); } + + template std::complex pmul(const std::complex& x, const std::complex& y) const + { return std::complex(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); } +}; + +template<> struct ei_conj_helper +{ + template std::complex + pmadd(const std::complex& x, const std::complex& y, const std::complex& c) const + { return c + pmul(x,y); } + + template std::complex pmul(const std::complex& x, const std::complex& y) const + { return std::complex(ei_real(x)*ei_real(y) - ei_imag(x)*ei_imag(y), - ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); } +}; + +// lightweight helper class to access matrix coefficients +template +class ei_blas_data_mapper +{ + public: + ei_blas_data_mapper(Scalar* data, int stride) : m_data(data), m_stride(stride) {} + EIGEN_STRONG_INLINE Scalar& operator()(int i, int j) + { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; } + protected: + Scalar* EIGEN_RESTRICT m_data; + int m_stride; +}; + +// lightweight helper class to access matrix coefficients (const version) +template +class ei_const_blas_data_mapper +{ + public: + ei_const_blas_data_mapper(const Scalar* data, int stride) : m_data(data), m_stride(stride) {} + EIGEN_STRONG_INLINE const Scalar& operator()(int i, int j) const + { return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; } + protected: + const Scalar* EIGEN_RESTRICT m_data; + int m_stride; +}; + +// +// template +// struct ei_L2_block_traits { +// enum {width = 8 * ei_meta_sqrt::ret }; +// }; + +// Defines various constant controlling level 3 blocking +template +struct ei_product_blocking_traits +{ + typedef typename ei_packet_traits::type PacketType; + enum { + PacketSize = sizeof(PacketType)/sizeof(Scalar), + #if (defined __i386__) + HalfRegisterCount = 4, + #else + HalfRegisterCount = 8, + #endif + + // register block size along the N direction + nr = HalfRegisterCount/2, + + // register block size along the M direction + mr = 2 * PacketSize, + + // max cache block size along the K direction + Max_kc = 8 * ei_meta_sqrt::ret, + + // max cache block size along the M direction + Max_mc = 2*Max_kc + }; +}; + +#endif // EIGEN_BLASUTIL_H