| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008-2016 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_GENERAL_MATRIX_VECTOR_H |
| #define EIGEN_GENERAL_MATRIX_VECTOR_H |
| |
| namespace Eigen { |
| |
| namespace internal { |
| |
| enum GEMVPacketSizeType { |
| GEMVPacketFull = 0, |
| GEMVPacketHalf, |
| GEMVPacketQuarter |
| }; |
| |
| template <int N, typename T1, typename T2, typename T3> |
| struct gemv_packet_cond { typedef T3 type; }; |
| |
| template <typename T1, typename T2, typename T3> |
| struct gemv_packet_cond<GEMVPacketFull, T1, T2, T3> { typedef T1 type; }; |
| |
| template <typename T1, typename T2, typename T3> |
| struct gemv_packet_cond<GEMVPacketHalf, T1, T2, T3> { typedef T2 type; }; |
| |
| template<typename LhsScalar, typename RhsScalar, int _PacketSize=GEMVPacketFull> |
| class gemv_traits |
| { |
| typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar; |
| |
| #define PACKET_DECL_COND_PREFIX(prefix, name, packet_size) \ |
| typedef typename gemv_packet_cond<packet_size, \ |
| typename packet_traits<name ## Scalar>::type, \ |
| typename packet_traits<name ## Scalar>::half, \ |
| typename unpacket_traits<typename packet_traits<name ## Scalar>::half>::half>::type \ |
| prefix ## name ## Packet |
| |
| PACKET_DECL_COND_PREFIX(_, Lhs, _PacketSize); |
| PACKET_DECL_COND_PREFIX(_, Rhs, _PacketSize); |
| PACKET_DECL_COND_PREFIX(_, Res, _PacketSize); |
| #undef PACKET_DECL_COND_PREFIX |
| |
| public: |
| enum { |
| Vectorizable = unpacket_traits<_LhsPacket>::vectorizable && |
| unpacket_traits<_RhsPacket>::vectorizable && |
| int(unpacket_traits<_LhsPacket>::size)==int(unpacket_traits<_RhsPacket>::size), |
| LhsPacketSize = Vectorizable ? unpacket_traits<_LhsPacket>::size : 1, |
| RhsPacketSize = Vectorizable ? unpacket_traits<_RhsPacket>::size : 1, |
| ResPacketSize = Vectorizable ? unpacket_traits<_ResPacket>::size : 1 |
| }; |
| |
| typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket; |
| typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket; |
| typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket; |
| }; |
| |
| |
| /* Optimized col-major matrix * vector product: |
| * This algorithm processes the matrix per vertical panels, |
| * which are then processed horizontaly per chunck of 8*PacketSize x 1 vertical segments. |
| * |
| * Mixing type logic: C += alpha * A * B |
| * | A | B |alpha| comments |
| * |real |cplx |cplx | no vectorization |
| * |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization |
| * |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp |
| * |cplx |real |real | optimal case, vectorization possible via real-cplx mul |
| * |
| * The same reasoning apply for the transposed case. |
| */ |
| template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version> |
| struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version> |
| { |
| typedef gemv_traits<LhsScalar,RhsScalar> Traits; |
| typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketHalf> HalfTraits; |
| typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketQuarter> QuarterTraits; |
| |
| typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar; |
| |
| typedef typename Traits::LhsPacket LhsPacket; |
| typedef typename Traits::RhsPacket RhsPacket; |
| typedef typename Traits::ResPacket ResPacket; |
| |
| typedef typename HalfTraits::LhsPacket LhsPacketHalf; |
| typedef typename HalfTraits::RhsPacket RhsPacketHalf; |
| typedef typename HalfTraits::ResPacket ResPacketHalf; |
| |
| typedef typename QuarterTraits::LhsPacket LhsPacketQuarter; |
| typedef typename QuarterTraits::RhsPacket RhsPacketQuarter; |
| typedef typename QuarterTraits::ResPacket ResPacketQuarter; |
| |
| EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run( |
| Index rows, Index cols, |
| const LhsMapper& lhs, |
| const RhsMapper& rhs, |
| ResScalar* res, Index resIncr, |
| RhsScalar alpha); |
| }; |
| |
| template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version> |
| EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run( |
| Index rows, Index cols, |
| const LhsMapper& alhs, |
| const RhsMapper& rhs, |
| ResScalar* res, Index resIncr, |
| RhsScalar alpha) |
| { |
| EIGEN_UNUSED_VARIABLE(resIncr); |
| eigen_internal_assert(resIncr==1); |
| |
| // The following copy tells the compiler that lhs's attributes are not modified outside this function |
| // This helps GCC to generate propoer code. |
| LhsMapper lhs(alhs); |
| |
| conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj; |
| conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj; |
| conj_helper<LhsPacketHalf,RhsPacketHalf,ConjugateLhs,ConjugateRhs> pcj_half; |
| conj_helper<LhsPacketQuarter,RhsPacketQuarter,ConjugateLhs,ConjugateRhs> pcj_quarter; |
| |
| const Index lhsStride = lhs.stride(); |
| // TODO: for padded aligned inputs, we could enable aligned reads |
| enum { LhsAlignment = Unaligned, |
| ResPacketSize = Traits::ResPacketSize, |
| ResPacketSizeHalf = HalfTraits::ResPacketSize, |
| ResPacketSizeQuarter = QuarterTraits::ResPacketSize, |
| LhsPacketSize = Traits::LhsPacketSize, |
| HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize, |
| HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf |
| }; |
| |
| const Index n8 = rows-8*ResPacketSize+1; |
| const Index n4 = rows-4*ResPacketSize+1; |
| const Index n3 = rows-3*ResPacketSize+1; |
| const Index n2 = rows-2*ResPacketSize+1; |
| const Index n1 = rows-1*ResPacketSize+1; |
| const Index n_half = rows-1*ResPacketSizeHalf+1; |
| const Index n_quarter = rows-1*ResPacketSizeQuarter+1; |
| |
| // TODO: improve the following heuristic: |
| const Index block_cols = cols<128 ? cols : (lhsStride*sizeof(LhsScalar)<32000?16:4); |
| ResPacket palpha = pset1<ResPacket>(alpha); |
| ResPacketHalf palpha_half = pset1<ResPacketHalf>(alpha); |
| ResPacketQuarter palpha_quarter = pset1<ResPacketQuarter>(alpha); |
| |
| for(Index j2=0; j2<cols; j2+=block_cols) |
| { |
| Index jend = numext::mini(j2+block_cols,cols); |
| Index i=0; |
| for(; i<n8; i+=ResPacketSize*8) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)), |
| c2 = pset1<ResPacket>(ResScalar(0)), |
| c3 = pset1<ResPacket>(ResScalar(0)), |
| c4 = pset1<ResPacket>(ResScalar(0)), |
| c5 = pset1<ResPacket>(ResScalar(0)), |
| c6 = pset1<ResPacket>(ResScalar(0)), |
| c7 = pset1<ResPacket>(ResScalar(0)); |
| |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacket b0 = pset1<RhsPacket>(rhs(j,0)); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1); |
| c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2); |
| c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*3,j),b0,c3); |
| c4 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*4,j),b0,c4); |
| c5 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*5,j),b0,c5); |
| c6 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*6,j),b0,c6); |
| c7 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*7,j),b0,c7); |
| } |
| pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0))); |
| pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1))); |
| pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2))); |
| pstoreu(res+i+ResPacketSize*3, pmadd(c3,palpha,ploadu<ResPacket>(res+i+ResPacketSize*3))); |
| pstoreu(res+i+ResPacketSize*4, pmadd(c4,palpha,ploadu<ResPacket>(res+i+ResPacketSize*4))); |
| pstoreu(res+i+ResPacketSize*5, pmadd(c5,palpha,ploadu<ResPacket>(res+i+ResPacketSize*5))); |
| pstoreu(res+i+ResPacketSize*6, pmadd(c6,palpha,ploadu<ResPacket>(res+i+ResPacketSize*6))); |
| pstoreu(res+i+ResPacketSize*7, pmadd(c7,palpha,ploadu<ResPacket>(res+i+ResPacketSize*7))); |
| } |
| if(i<n4) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)), |
| c2 = pset1<ResPacket>(ResScalar(0)), |
| c3 = pset1<ResPacket>(ResScalar(0)); |
| |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacket b0 = pset1<RhsPacket>(rhs(j,0)); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1); |
| c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2); |
| c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*3,j),b0,c3); |
| } |
| pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0))); |
| pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1))); |
| pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2))); |
| pstoreu(res+i+ResPacketSize*3, pmadd(c3,palpha,ploadu<ResPacket>(res+i+ResPacketSize*3))); |
| |
| i+=ResPacketSize*4; |
| } |
| if(i<n3) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)), |
| c2 = pset1<ResPacket>(ResScalar(0)); |
| |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacket b0 = pset1<RhsPacket>(rhs(j,0)); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1); |
| c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2); |
| } |
| pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0))); |
| pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1))); |
| pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2))); |
| |
| i+=ResPacketSize*3; |
| } |
| if(i<n2) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)); |
| |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacket b0 = pset1<RhsPacket>(rhs(j,0)); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1); |
| } |
| pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0))); |
| pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1))); |
| i+=ResPacketSize*2; |
| } |
| if(i<n1) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)); |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacket b0 = pset1<RhsPacket>(rhs(j,0)); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0); |
| } |
| pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0))); |
| i+=ResPacketSize; |
| } |
| if(HasHalf && i<n_half) |
| { |
| ResPacketHalf c0 = pset1<ResPacketHalf>(ResScalar(0)); |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacketHalf b0 = pset1<RhsPacketHalf>(rhs(j,0)); |
| c0 = pcj_half.pmadd(lhs.template load<LhsPacketHalf,LhsAlignment>(i+0,j),b0,c0); |
| } |
| pstoreu(res+i+ResPacketSizeHalf*0, pmadd(c0,palpha_half,ploadu<ResPacketHalf>(res+i+ResPacketSizeHalf*0))); |
| i+=ResPacketSizeHalf; |
| } |
| if(HasQuarter && i<n_quarter) |
| { |
| ResPacketQuarter c0 = pset1<ResPacketQuarter>(ResScalar(0)); |
| for(Index j=j2; j<jend; j+=1) |
| { |
| RhsPacketQuarter b0 = pset1<RhsPacketQuarter>(rhs(j,0)); |
| c0 = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter,LhsAlignment>(i+0,j),b0,c0); |
| } |
| pstoreu(res+i+ResPacketSizeQuarter*0, pmadd(c0,palpha_quarter,ploadu<ResPacketQuarter>(res+i+ResPacketSizeQuarter*0))); |
| i+=ResPacketSizeQuarter; |
| } |
| for(;i<rows;++i) |
| { |
| ResScalar c0(0); |
| for(Index j=j2; j<jend; j+=1) |
| c0 += cj.pmul(lhs(i,j), rhs(j,0)); |
| res[i] += alpha*c0; |
| } |
| } |
| } |
| |
| /* Optimized row-major matrix * vector product: |
| * This algorithm processes 4 rows at once that allows to both reduce |
| * the number of load/stores of the result by a factor 4 and to reduce |
| * the instruction dependency. Moreover, we know that all bands have the |
| * same alignment pattern. |
| * |
| * Mixing type logic: |
| * - alpha is always a complex (or converted to a complex) |
| * - no vectorization |
| */ |
| template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version> |
| struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version> |
| { |
| typedef gemv_traits<LhsScalar,RhsScalar> Traits; |
| typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketHalf> HalfTraits; |
| typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketQuarter> QuarterTraits; |
| |
| typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar; |
| |
| typedef typename Traits::LhsPacket LhsPacket; |
| typedef typename Traits::RhsPacket RhsPacket; |
| typedef typename Traits::ResPacket ResPacket; |
| |
| typedef typename HalfTraits::LhsPacket LhsPacketHalf; |
| typedef typename HalfTraits::RhsPacket RhsPacketHalf; |
| typedef typename HalfTraits::ResPacket ResPacketHalf; |
| |
| typedef typename QuarterTraits::LhsPacket LhsPacketQuarter; |
| typedef typename QuarterTraits::RhsPacket RhsPacketQuarter; |
| typedef typename QuarterTraits::ResPacket ResPacketQuarter; |
| |
| EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run( |
| Index rows, Index cols, |
| const LhsMapper& lhs, |
| const RhsMapper& rhs, |
| ResScalar* res, Index resIncr, |
| ResScalar alpha); |
| }; |
| |
| template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version> |
| EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run( |
| Index rows, Index cols, |
| const LhsMapper& alhs, |
| const RhsMapper& rhs, |
| ResScalar* res, Index resIncr, |
| ResScalar alpha) |
| { |
| // The following copy tells the compiler that lhs's attributes are not modified outside this function |
| // This helps GCC to generate propoer code. |
| LhsMapper lhs(alhs); |
| |
| eigen_internal_assert(rhs.stride()==1); |
| conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj; |
| conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj; |
| conj_helper<LhsPacketHalf,RhsPacketHalf,ConjugateLhs,ConjugateRhs> pcj_half; |
| conj_helper<LhsPacketQuarter,RhsPacketQuarter,ConjugateLhs,ConjugateRhs> pcj_quarter; |
| |
| // TODO: fine tune the following heuristic. The rationale is that if the matrix is very large, |
| // processing 8 rows at once might be counter productive wrt cache. |
| const Index n8 = lhs.stride()*sizeof(LhsScalar)>32000 ? 0 : rows-7; |
| const Index n4 = rows-3; |
| const Index n2 = rows-1; |
| |
| // TODO: for padded aligned inputs, we could enable aligned reads |
| enum { LhsAlignment = Unaligned, |
| ResPacketSize = Traits::ResPacketSize, |
| ResPacketSizeHalf = HalfTraits::ResPacketSize, |
| ResPacketSizeQuarter = QuarterTraits::ResPacketSize, |
| LhsPacketSize = Traits::LhsPacketSize, |
| LhsPacketSizeHalf = HalfTraits::LhsPacketSize, |
| LhsPacketSizeQuarter = QuarterTraits::LhsPacketSize, |
| HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize, |
| HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf |
| }; |
| |
| Index i=0; |
| for(; i<n8; i+=8) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)), |
| c2 = pset1<ResPacket>(ResScalar(0)), |
| c3 = pset1<ResPacket>(ResScalar(0)), |
| c4 = pset1<ResPacket>(ResScalar(0)), |
| c5 = pset1<ResPacket>(ResScalar(0)), |
| c6 = pset1<ResPacket>(ResScalar(0)), |
| c7 = pset1<ResPacket>(ResScalar(0)); |
| |
| Index j=0; |
| for(; j+LhsPacketSize<=cols; j+=LhsPacketSize) |
| { |
| RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0); |
| |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1); |
| c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+2,j),b0,c2); |
| c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+3,j),b0,c3); |
| c4 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+4,j),b0,c4); |
| c5 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+5,j),b0,c5); |
| c6 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+6,j),b0,c6); |
| c7 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+7,j),b0,c7); |
| } |
| ResScalar cc0 = predux(c0); |
| ResScalar cc1 = predux(c1); |
| ResScalar cc2 = predux(c2); |
| ResScalar cc3 = predux(c3); |
| ResScalar cc4 = predux(c4); |
| ResScalar cc5 = predux(c5); |
| ResScalar cc6 = predux(c6); |
| ResScalar cc7 = predux(c7); |
| for(; j<cols; ++j) |
| { |
| RhsScalar b0 = rhs(j,0); |
| |
| cc0 += cj.pmul(lhs(i+0,j), b0); |
| cc1 += cj.pmul(lhs(i+1,j), b0); |
| cc2 += cj.pmul(lhs(i+2,j), b0); |
| cc3 += cj.pmul(lhs(i+3,j), b0); |
| cc4 += cj.pmul(lhs(i+4,j), b0); |
| cc5 += cj.pmul(lhs(i+5,j), b0); |
| cc6 += cj.pmul(lhs(i+6,j), b0); |
| cc7 += cj.pmul(lhs(i+7,j), b0); |
| } |
| res[(i+0)*resIncr] += alpha*cc0; |
| res[(i+1)*resIncr] += alpha*cc1; |
| res[(i+2)*resIncr] += alpha*cc2; |
| res[(i+3)*resIncr] += alpha*cc3; |
| res[(i+4)*resIncr] += alpha*cc4; |
| res[(i+5)*resIncr] += alpha*cc5; |
| res[(i+6)*resIncr] += alpha*cc6; |
| res[(i+7)*resIncr] += alpha*cc7; |
| } |
| for(; i<n4; i+=4) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)), |
| c2 = pset1<ResPacket>(ResScalar(0)), |
| c3 = pset1<ResPacket>(ResScalar(0)); |
| |
| Index j=0; |
| for(; j+LhsPacketSize<=cols; j+=LhsPacketSize) |
| { |
| RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0); |
| |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1); |
| c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+2,j),b0,c2); |
| c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+3,j),b0,c3); |
| } |
| ResScalar cc0 = predux(c0); |
| ResScalar cc1 = predux(c1); |
| ResScalar cc2 = predux(c2); |
| ResScalar cc3 = predux(c3); |
| for(; j<cols; ++j) |
| { |
| RhsScalar b0 = rhs(j,0); |
| |
| cc0 += cj.pmul(lhs(i+0,j), b0); |
| cc1 += cj.pmul(lhs(i+1,j), b0); |
| cc2 += cj.pmul(lhs(i+2,j), b0); |
| cc3 += cj.pmul(lhs(i+3,j), b0); |
| } |
| res[(i+0)*resIncr] += alpha*cc0; |
| res[(i+1)*resIncr] += alpha*cc1; |
| res[(i+2)*resIncr] += alpha*cc2; |
| res[(i+3)*resIncr] += alpha*cc3; |
| } |
| for(; i<n2; i+=2) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)), |
| c1 = pset1<ResPacket>(ResScalar(0)); |
| |
| Index j=0; |
| for(; j+LhsPacketSize<=cols; j+=LhsPacketSize) |
| { |
| RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0); |
| |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0); |
| c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1); |
| } |
| ResScalar cc0 = predux(c0); |
| ResScalar cc1 = predux(c1); |
| for(; j<cols; ++j) |
| { |
| RhsScalar b0 = rhs(j,0); |
| |
| cc0 += cj.pmul(lhs(i+0,j), b0); |
| cc1 += cj.pmul(lhs(i+1,j), b0); |
| } |
| res[(i+0)*resIncr] += alpha*cc0; |
| res[(i+1)*resIncr] += alpha*cc1; |
| } |
| for(; i<rows; ++i) |
| { |
| ResPacket c0 = pset1<ResPacket>(ResScalar(0)); |
| ResPacketHalf c0_h = pset1<ResPacketHalf>(ResScalar(0)); |
| ResPacketQuarter c0_q = pset1<ResPacketQuarter>(ResScalar(0)); |
| Index j=0; |
| for(; j+LhsPacketSize<=cols; j+=LhsPacketSize) |
| { |
| RhsPacket b0 = rhs.template load<RhsPacket,Unaligned>(j,0); |
| c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i,j),b0,c0); |
| } |
| ResScalar cc0 = predux(c0); |
| if (HasHalf) { |
| for(; j+LhsPacketSizeHalf<=cols; j+=LhsPacketSizeHalf) |
| { |
| RhsPacketHalf b0 = rhs.template load<RhsPacketHalf,Unaligned>(j,0); |
| c0_h = pcj_half.pmadd(lhs.template load<LhsPacketHalf,LhsAlignment>(i,j),b0,c0_h); |
| } |
| cc0 += predux(c0_h); |
| } |
| if (HasQuarter) { |
| for(; j+LhsPacketSizeQuarter<=cols; j+=LhsPacketSizeQuarter) |
| { |
| RhsPacketQuarter b0 = rhs.template load<RhsPacketQuarter,Unaligned>(j,0); |
| c0_q = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter,LhsAlignment>(i,j),b0,c0_q); |
| } |
| cc0 += predux(c0_q); |
| } |
| for(; j<cols; ++j) |
| { |
| cc0 += cj.pmul(lhs(i,j), rhs(j,0)); |
| } |
| res[i*resIncr] += alpha*cc0; |
| } |
| } |
| |
| } // end namespace internal |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_GENERAL_MATRIX_VECTOR_H |
