Quelle GeneralMatrixVector.h Sprache: C

// 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

quality80%

¤ Dauer der Verarbeitung: 0.16 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.