GeneralBlockPanelKernel.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_GENERAL_BLOCK_PANEL_H
26 #define EIGEN_GENERAL_BLOCK_PANEL_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
32 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs=false, bool _ConjRhs=false>
33 class gebp_traits;
34 
35 
37 inline std::ptrdiff_t manage_caching_sizes_helper(std::ptrdiff_t a, std::ptrdiff_t b)
38 {
39  return a<=0 ? b : a;
40 }
41 
43 inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdiff_t* l2=0)
44 {
45  static std::ptrdiff_t m_l1CacheSize = 0;
46  static std::ptrdiff_t m_l2CacheSize = 0;
47  if(m_l2CacheSize==0)
48  {
49  m_l1CacheSize = manage_caching_sizes_helper(queryL1CacheSize(),8 * 1024);
50  m_l2CacheSize = manage_caching_sizes_helper(queryTopLevelCacheSize(),1*1024*1024);
51  }
52 
53  if(action==SetAction)
54  {
55  // set the cpu cache size and cache all block sizes from a global cache size in byte
56  eigen_internal_assert(l1!=0 && l2!=0);
57  m_l1CacheSize = *l1;
58  m_l2CacheSize = *l2;
59  }
60  else if(action==GetAction)
61  {
62  eigen_internal_assert(l1!=0 && l2!=0);
63  *l1 = m_l1CacheSize;
64  *l2 = m_l2CacheSize;
65  }
66  else
67  {
68  eigen_internal_assert(false);
69  }
70 }
71 
87 template<typename LhsScalar, typename RhsScalar, int KcFactor>
88 void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrdiff_t& n)
89 {
91  // Explanations:
92  // Let's recall the product algorithms form kc x nc horizontal panels B' on the rhs and
93  // mc x kc blocks A' on the lhs. A' has to fit into L2 cache. Moreover, B' is processed
94  // per kc x nr vertical small panels where nr is the blocking size along the n dimension
95  // at the register level. For vectorization purpose, these small vertical panels are unpacked,
96  // e.g., each coefficient is replicated to fit a packet. This small vertical panel has to
97  // stay in L1 cache.
98  std::ptrdiff_t l1, l2;
99 
100  typedef gebp_traits<LhsScalar,RhsScalar> Traits;
101  enum {
102  kdiv = KcFactor * 2 * Traits::nr
103  * Traits::RhsProgress * sizeof(RhsScalar),
104  mr = gebp_traits<LhsScalar,RhsScalar>::mr,
105  mr_mask = (0xffffffff/mr)*mr
106  };
107 
108  manage_caching_sizes(GetAction, &l1, &l2);
109  k = std::min<std::ptrdiff_t>(k, l1/kdiv);
110  std::ptrdiff_t _m = k>0 ? l2/(4 * sizeof(LhsScalar) * k) : 0;
111  if(_m<m) m = _m & mr_mask;
112 }
113 
114 template<typename LhsScalar, typename RhsScalar>
115 inline void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrdiff_t& n)
116 {
117  computeProductBlockingSizes<LhsScalar,RhsScalar,1>(k, m, n);
118 }
119 
120 #ifdef EIGEN_HAS_FUSE_CJMADD
121  #define MADD(CJ,A,B,C,T) C = CJ.pmadd(A,B,C);
122 #else
123 
124  // FIXME (a bit overkill maybe ?)
125 
126  template<typename CJ, typename A, typename B, typename C, typename T> struct gebp_madd_selector {
127  EIGEN_ALWAYS_INLINE static void run(const CJ& cj, A& a, B& b, C& c, T& /*t*/)
128  {
129  c = cj.pmadd(a,b,c);
130  }
131  };
132 
133  template<typename CJ, typename T> struct gebp_madd_selector<CJ,T,T,T,T> {
134  EIGEN_ALWAYS_INLINE static void run(const CJ& cj, T& a, T& b, T& c, T& t)
135  {
136  t = b; t = cj.pmul(a,t); c = padd(c,t);
137  }
138  };
139 
140  template<typename CJ, typename A, typename B, typename C, typename T>
141  EIGEN_STRONG_INLINE void gebp_madd(const CJ& cj, A& a, B& b, C& c, T& t)
142  {
143  gebp_madd_selector<CJ,A,B,C,T>::run(cj,a,b,c,t);
144  }
145 
146  #define MADD(CJ,A,B,C,T) gebp_madd(CJ,A,B,C,T);
147 // #define MADD(CJ,A,B,C,T) T = B; T = CJ.pmul(A,T); C = padd(C,T);
148 #endif
149 
150 /* Vectorization logic
151  * real*real: unpack rhs to constant packets, ...
152  *
153  * cd*cd : unpack rhs to (b_r,b_r), (b_i,b_i), mul to get (a_r b_r,a_i b_r) (a_r b_i,a_i b_i),
154  * storing each res packet into two packets (2x2),
155  * at the end combine them: swap the second and addsub them
156  * cf*cf : same but with 2x4 blocks
157  * cplx*real : unpack rhs to constant packets, ...
158  * real*cplx : load lhs as (a0,a0,a1,a1), and mul as usual
159  */
160 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs, bool _ConjRhs>
161 class gebp_traits
162 {
163 public:
164  typedef _LhsScalar LhsScalar;
165  typedef _RhsScalar RhsScalar;
166  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
167 
168  enum {
169  ConjLhs = _ConjLhs,
170  ConjRhs = _ConjRhs,
171  Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable,
172  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
173  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
174  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
175 
176  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
177 
178  // register block size along the N direction (must be either 2 or 4)
179  nr = NumberOfRegisters/4,
180 
181  // register block size along the M direction (currently, this one cannot be modified)
182  mr = 2 * LhsPacketSize,
183 
184  WorkSpaceFactor = nr * RhsPacketSize,
185 
186  LhsProgress = LhsPacketSize,
187  RhsProgress = RhsPacketSize
188  };
189 
190  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
191  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
192  typedef typename packet_traits<ResScalar>::type _ResPacket;
193 
194  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
195  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
196  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
197 
198  typedef ResPacket AccPacket;
199 
200  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
201  {
202  p = pset1<ResPacket>(ResScalar(0));
203  }
204 
205  EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
206  {
207  for(DenseIndex k=0; k<n; k++)
208  pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
209  }
210 
211  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
212  {
213  dest = pload<RhsPacket>(b);
214  }
215 
216  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
217  {
218  dest = pload<LhsPacket>(a);
219  }
220 
221  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, AccPacket& c, AccPacket& tmp) const
222  {
223  tmp = b; tmp = pmul(a,tmp); c = padd(c,tmp);
224  }
225 
226  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
227  {
228  r = pmadd(c,alpha,r);
229  }
230 
231 protected:
232 // conj_helper<LhsScalar,RhsScalar,ConjLhs,ConjRhs> cj;
233 // conj_helper<LhsPacket,RhsPacket,ConjLhs,ConjRhs> pcj;
234 };
235 
236 template<typename RealScalar, bool _ConjLhs>
237 class gebp_traits<std::complex<RealScalar>, RealScalar, _ConjLhs, false>
238 {
239 public:
240  typedef std::complex<RealScalar> LhsScalar;
241  typedef RealScalar RhsScalar;
242  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
243 
244  enum {
245  ConjLhs = _ConjLhs,
246  ConjRhs = false,
247  Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable,
248  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
249  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
250  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
251 
252  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
253  nr = NumberOfRegisters/4,
254  mr = 2 * LhsPacketSize,
255  WorkSpaceFactor = nr*RhsPacketSize,
256 
257  LhsProgress = LhsPacketSize,
258  RhsProgress = RhsPacketSize
259  };
260 
261  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
262  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
263  typedef typename packet_traits<ResScalar>::type _ResPacket;
264 
265  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
266  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
267  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
268 
269  typedef ResPacket AccPacket;
270 
271  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
272  {
273  p = pset1<ResPacket>(ResScalar(0));
274  }
275 
276  EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
277  {
278  for(DenseIndex k=0; k<n; k++)
279  pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
280  }
281 
282  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
283  {
284  dest = pload<RhsPacket>(b);
285  }
286 
287  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
288  {
289  dest = pload<LhsPacket>(a);
290  }
291 
292  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp) const
293  {
294  madd_impl(a, b, c, tmp, typename conditional<Vectorizable,true_type,false_type>::type());
295  }
296 
297  EIGEN_STRONG_INLINE void madd_impl(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp, const true_type&) const
298  {
299  tmp = b; tmp = pmul(a.v,tmp); c.v = padd(c.v,tmp);
300  }
301 
302  EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/, const false_type&) const
303  {
304  c += a * b;
305  }
306 
307  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
308  {
309  r = cj.pmadd(c,alpha,r);
310  }
311 
312 protected:
313  conj_helper<ResPacket,ResPacket,ConjLhs,false> cj;
314 };
315 
316 template<typename RealScalar, bool _ConjLhs, bool _ConjRhs>
317 class gebp_traits<std::complex<RealScalar>, std::complex<RealScalar>, _ConjLhs, _ConjRhs >
318 {
319 public:
320  typedef std::complex<RealScalar> Scalar;
321  typedef std::complex<RealScalar> LhsScalar;
322  typedef std::complex<RealScalar> RhsScalar;
323  typedef std::complex<RealScalar> ResScalar;
324 
325  enum {
326  ConjLhs = _ConjLhs,
327  ConjRhs = _ConjRhs,
328  Vectorizable = packet_traits<RealScalar>::Vectorizable
329  && packet_traits<Scalar>::Vectorizable,
330  RealPacketSize = Vectorizable ? packet_traits<RealScalar>::size : 1,
331  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
332 
333  nr = 2,
334  mr = 2 * ResPacketSize,
335  WorkSpaceFactor = Vectorizable ? 2*nr*RealPacketSize : nr,
336 
337  LhsProgress = ResPacketSize,
338  RhsProgress = Vectorizable ? 2*ResPacketSize : 1
339  };
340 
341  typedef typename packet_traits<RealScalar>::type RealPacket;
342  typedef typename packet_traits<Scalar>::type ScalarPacket;
343  struct DoublePacket
344  {
345  RealPacket first;
346  RealPacket second;
347  };
348 
349  typedef typename conditional<Vectorizable,RealPacket, Scalar>::type LhsPacket;
350  typedef typename conditional<Vectorizable,DoublePacket,Scalar>::type RhsPacket;
351  typedef typename conditional<Vectorizable,ScalarPacket,Scalar>::type ResPacket;
352  typedef typename conditional<Vectorizable,DoublePacket,Scalar>::type AccPacket;
353 
354  EIGEN_STRONG_INLINE void initAcc(Scalar& p) { p = Scalar(0); }
355 
356  EIGEN_STRONG_INLINE void initAcc(DoublePacket& p)
357  {
358  p.first = pset1<RealPacket>(RealScalar(0));
359  p.second = pset1<RealPacket>(RealScalar(0));
360  }
361 
362  /* Unpack the rhs coeff such that each complex coefficient is spread into
363  * two packects containing respectively the real and imaginary coefficient
364  * duplicated as many time as needed: (x+iy) => [x, ..., x] [y, ..., y]
365  */
366  EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const Scalar* rhs, Scalar* b)
367  {
368  for(DenseIndex k=0; k<n; k++)
369  {
370  if(Vectorizable)
371  {
372  pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+0], real(rhs[k]));
373  pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+ResPacketSize], imag(rhs[k]));
374  }
375  else
376  b[k] = rhs[k];
377  }
378  }
379 
380  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, ResPacket& dest) const { dest = *b; }
381 
382  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, DoublePacket& dest) const
383  {
384  dest.first = pload<RealPacket>((const RealScalar*)b);
385  dest.second = pload<RealPacket>((const RealScalar*)(b+ResPacketSize));
386  }
387 
388  // nothing special here
389  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
390  {
391  dest = pload<LhsPacket>((const typename unpacket_traits<LhsPacket>::type*)(a));
392  }
393 
394  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, DoublePacket& c, RhsPacket& /*tmp*/) const
395  {
396  c.first = padd(pmul(a,b.first), c.first);
397  c.second = padd(pmul(a,b.second),c.second);
398  }
399 
400  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, ResPacket& c, RhsPacket& /*tmp*/) const
401  {
402  c = cj.pmadd(a,b,c);
403  }
404 
405  EIGEN_STRONG_INLINE void acc(const Scalar& c, const Scalar& alpha, Scalar& r) const { r += alpha * c; }
406 
407  EIGEN_STRONG_INLINE void acc(const DoublePacket& c, const ResPacket& alpha, ResPacket& r) const
408  {
409  // assemble c
410  ResPacket tmp;
411  if((!ConjLhs)&&(!ConjRhs))
412  {
413  tmp = pcplxflip(pconj(ResPacket(c.second)));
414  tmp = padd(ResPacket(c.first),tmp);
415  }
416  else if((!ConjLhs)&&(ConjRhs))
417  {
418  tmp = pconj(pcplxflip(ResPacket(c.second)));
419  tmp = padd(ResPacket(c.first),tmp);
420  }
421  else if((ConjLhs)&&(!ConjRhs))
422  {
423  tmp = pcplxflip(ResPacket(c.second));
424  tmp = padd(pconj(ResPacket(c.first)),tmp);
425  }
426  else if((ConjLhs)&&(ConjRhs))
427  {
428  tmp = pcplxflip(ResPacket(c.second));
429  tmp = psub(pconj(ResPacket(c.first)),tmp);
430  }
431 
432  r = pmadd(tmp,alpha,r);
433  }
434 
435 protected:
436  conj_helper<LhsScalar,RhsScalar,ConjLhs,ConjRhs> cj;
437 };
438 
439 template<typename RealScalar, bool _ConjRhs>
440 class gebp_traits<RealScalar, std::complex<RealScalar>, false, _ConjRhs >
441 {
442 public:
443  typedef std::complex<RealScalar> Scalar;
444  typedef RealScalar LhsScalar;
445  typedef Scalar RhsScalar;
446  typedef Scalar ResScalar;
447 
448  enum {
449  ConjLhs = false,
450  ConjRhs = _ConjRhs,
451  Vectorizable = packet_traits<RealScalar>::Vectorizable
452  && packet_traits<Scalar>::Vectorizable,
453  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
454  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
455  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
456 
457  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
458  nr = 4,
459  mr = 2*ResPacketSize,
460  WorkSpaceFactor = nr*RhsPacketSize,
461 
462  LhsProgress = ResPacketSize,
463  RhsProgress = ResPacketSize
464  };
465 
466  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
467  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
468  typedef typename packet_traits<ResScalar>::type _ResPacket;
469 
470  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
471  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
472  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
473 
474  typedef ResPacket AccPacket;
475 
476  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
477  {
478  p = pset1<ResPacket>(ResScalar(0));
479  }
480 
481  EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
482  {
483  for(DenseIndex k=0; k<n; k++)
484  pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
485  }
486 
487  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
488  {
489  dest = pload<RhsPacket>(b);
490  }
491 
492  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
493  {
494  dest = ploaddup<LhsPacket>(a);
495  }
496 
497  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp) const
498  {
499  madd_impl(a, b, c, tmp, typename conditional<Vectorizable,true_type,false_type>::type());
500  }
501 
502  EIGEN_STRONG_INLINE void madd_impl(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp, const true_type&) const
503  {
504  tmp = b; tmp.v = pmul(a,tmp.v); c = padd(c,tmp);
505  }
506 
507  EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/, const false_type&) const
508  {
509  c += a * b;
510  }
511 
512  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
513  {
514  r = cj.pmadd(alpha,c,r);
515  }
516 
517 protected:
518  conj_helper<ResPacket,ResPacket,false,ConjRhs> cj;
519 };
520 
521 /* optimized GEneral packed Block * packed Panel product kernel
522  *
523  * Mixing type logic: C += A * B
524  * | A | B | comments
525  * |real |cplx | no vectorization yet, would require to pack A with duplication
526  * |cplx |real | easy vectorization
527  */
528 template<typename LhsScalar, typename RhsScalar, typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
529 struct gebp_kernel
530 {
531  typedef gebp_traits<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> Traits;
532  typedef typename Traits::ResScalar ResScalar;
533  typedef typename Traits::LhsPacket LhsPacket;
534  typedef typename Traits::RhsPacket RhsPacket;
535  typedef typename Traits::ResPacket ResPacket;
536  typedef typename Traits::AccPacket AccPacket;
537 
538  enum {
539  Vectorizable = Traits::Vectorizable,
540  LhsProgress = Traits::LhsProgress,
541  RhsProgress = Traits::RhsProgress,
542  ResPacketSize = Traits::ResPacketSize
543  };
544 
546  void operator()(ResScalar* res, Index resStride, const LhsScalar* blockA, const RhsScalar* blockB, Index rows, Index depth, Index cols, ResScalar alpha,
547  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, RhsScalar* unpackedB = 0)
548  {
549  Traits traits;
550 
551  if(strideA==-1) strideA = depth;
552  if(strideB==-1) strideB = depth;
553  conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
554 // conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj;
555  Index packet_cols = (cols/nr) * nr;
556  const Index peeled_mc = (rows/mr)*mr;
557  // FIXME:
558  const Index peeled_mc2 = peeled_mc + (rows-peeled_mc >= LhsProgress ? LhsProgress : 0);
559  const Index peeled_kc = (depth/4)*4;
560 
561  if(unpackedB==0)
562  unpackedB = const_cast<RhsScalar*>(blockB - strideB * nr * RhsProgress);
563 
564  // loops on each micro vertical panel of rhs (depth x nr)
565  for(Index j2=0; j2<packet_cols; j2+=nr)
566  {
567  traits.unpackRhs(depth*nr,&blockB[j2*strideB+offsetB*nr],unpackedB);
568 
569  // loops on each largest micro horizontal panel of lhs (mr x depth)
570  // => we select a mr x nr micro block of res which is entirely
571  // stored into mr/packet_size x nr registers.
572  for(Index i=0; i<peeled_mc; i+=mr)
573  {
574  const LhsScalar* blA = &blockA[i*strideA+offsetA*mr];
575  prefetch(&blA[0]);
576 
577  // gets res block as register
578  AccPacket C0, C1, C2, C3, C4, C5, C6, C7;
579  traits.initAcc(C0);
580  traits.initAcc(C1);
581  if(nr==4) traits.initAcc(C2);
582  if(nr==4) traits.initAcc(C3);
583  traits.initAcc(C4);
584  traits.initAcc(C5);
585  if(nr==4) traits.initAcc(C6);
586  if(nr==4) traits.initAcc(C7);
587 
588  ResScalar* r0 = &res[(j2+0)*resStride + i];
589  ResScalar* r1 = r0 + resStride;
590  ResScalar* r2 = r1 + resStride;
591  ResScalar* r3 = r2 + resStride;
592 
593  prefetch(r0+16);
594  prefetch(r1+16);
595  prefetch(r2+16);
596  prefetch(r3+16);
597 
598  // performs "inner" product
599  // TODO let's check wether the folowing peeled loop could not be
600  // optimized via optimal prefetching from one loop to the other
601  const RhsScalar* blB = unpackedB;
602  for(Index k=0; k<peeled_kc; k+=4)
603  {
604  if(nr==2)
605  {
606  LhsPacket A0, A1;
607  RhsPacket B_0;
608  RhsPacket T0;
609 
610 EIGEN_ASM_COMMENT("mybegin2");
611  traits.loadLhs(&blA[0*LhsProgress], A0);
612  traits.loadLhs(&blA[1*LhsProgress], A1);
613  traits.loadRhs(&blB[0*RhsProgress], B_0);
614  traits.madd(A0,B_0,C0,T0);
615  traits.madd(A1,B_0,C4,B_0);
616  traits.loadRhs(&blB[1*RhsProgress], B_0);
617  traits.madd(A0,B_0,C1,T0);
618  traits.madd(A1,B_0,C5,B_0);
619 
620  traits.loadLhs(&blA[2*LhsProgress], A0);
621  traits.loadLhs(&blA[3*LhsProgress], A1);
622  traits.loadRhs(&blB[2*RhsProgress], B_0);
623  traits.madd(A0,B_0,C0,T0);
624  traits.madd(A1,B_0,C4,B_0);
625  traits.loadRhs(&blB[3*RhsProgress], B_0);
626  traits.madd(A0,B_0,C1,T0);
627  traits.madd(A1,B_0,C5,B_0);
628 
629  traits.loadLhs(&blA[4*LhsProgress], A0);
630  traits.loadLhs(&blA[5*LhsProgress], A1);
631  traits.loadRhs(&blB[4*RhsProgress], B_0);
632  traits.madd(A0,B_0,C0,T0);
633  traits.madd(A1,B_0,C4,B_0);
634  traits.loadRhs(&blB[5*RhsProgress], B_0);
635  traits.madd(A0,B_0,C1,T0);
636  traits.madd(A1,B_0,C5,B_0);
637 
638  traits.loadLhs(&blA[6*LhsProgress], A0);
639  traits.loadLhs(&blA[7*LhsProgress], A1);
640  traits.loadRhs(&blB[6*RhsProgress], B_0);
641  traits.madd(A0,B_0,C0,T0);
642  traits.madd(A1,B_0,C4,B_0);
643  traits.loadRhs(&blB[7*RhsProgress], B_0);
644  traits.madd(A0,B_0,C1,T0);
645  traits.madd(A1,B_0,C5,B_0);
646 EIGEN_ASM_COMMENT("myend");
647  }
648  else
649  {
650 EIGEN_ASM_COMMENT("mybegin4");
651  LhsPacket A0, A1;
652  RhsPacket B_0, B1, B2, B3;
653  RhsPacket T0;
654 
655  traits.loadLhs(&blA[0*LhsProgress], A0);
656  traits.loadLhs(&blA[1*LhsProgress], A1);
657  traits.loadRhs(&blB[0*RhsProgress], B_0);
658  traits.loadRhs(&blB[1*RhsProgress], B1);
659 
660  traits.madd(A0,B_0,C0,T0);
661  traits.loadRhs(&blB[2*RhsProgress], B2);
662  traits.madd(A1,B_0,C4,B_0);
663  traits.loadRhs(&blB[3*RhsProgress], B3);
664  traits.loadRhs(&blB[4*RhsProgress], B_0);
665  traits.madd(A0,B1,C1,T0);
666  traits.madd(A1,B1,C5,B1);
667  traits.loadRhs(&blB[5*RhsProgress], B1);
668  traits.madd(A0,B2,C2,T0);
669  traits.madd(A1,B2,C6,B2);
670  traits.loadRhs(&blB[6*RhsProgress], B2);
671  traits.madd(A0,B3,C3,T0);
672  traits.loadLhs(&blA[2*LhsProgress], A0);
673  traits.madd(A1,B3,C7,B3);
674  traits.loadLhs(&blA[3*LhsProgress], A1);
675  traits.loadRhs(&blB[7*RhsProgress], B3);
676  traits.madd(A0,B_0,C0,T0);
677  traits.madd(A1,B_0,C4,B_0);
678  traits.loadRhs(&blB[8*RhsProgress], B_0);
679  traits.madd(A0,B1,C1,T0);
680  traits.madd(A1,B1,C5,B1);
681  traits.loadRhs(&blB[9*RhsProgress], B1);
682  traits.madd(A0,B2,C2,T0);
683  traits.madd(A1,B2,C6,B2);
684  traits.loadRhs(&blB[10*RhsProgress], B2);
685  traits.madd(A0,B3,C3,T0);
686  traits.loadLhs(&blA[4*LhsProgress], A0);
687  traits.madd(A1,B3,C7,B3);
688  traits.loadLhs(&blA[5*LhsProgress], A1);
689  traits.loadRhs(&blB[11*RhsProgress], B3);
690 
691  traits.madd(A0,B_0,C0,T0);
692  traits.madd(A1,B_0,C4,B_0);
693  traits.loadRhs(&blB[12*RhsProgress], B_0);
694  traits.madd(A0,B1,C1,T0);
695  traits.madd(A1,B1,C5,B1);
696  traits.loadRhs(&blB[13*RhsProgress], B1);
697  traits.madd(A0,B2,C2,T0);
698  traits.madd(A1,B2,C6,B2);
699  traits.loadRhs(&blB[14*RhsProgress], B2);
700  traits.madd(A0,B3,C3,T0);
701  traits.loadLhs(&blA[6*LhsProgress], A0);
702  traits.madd(A1,B3,C7,B3);
703  traits.loadLhs(&blA[7*LhsProgress], A1);
704  traits.loadRhs(&blB[15*RhsProgress], B3);
705  traits.madd(A0,B_0,C0,T0);
706  traits.madd(A1,B_0,C4,B_0);
707  traits.madd(A0,B1,C1,T0);
708  traits.madd(A1,B1,C5,B1);
709  traits.madd(A0,B2,C2,T0);
710  traits.madd(A1,B2,C6,B2);
711  traits.madd(A0,B3,C3,T0);
712  traits.madd(A1,B3,C7,B3);
713  }
714 
715  blB += 4*nr*RhsProgress;
716  blA += 4*mr;
717  }
718  // process remaining peeled loop
719  for(Index k=peeled_kc; k<depth; k++)
720  {
721  if(nr==2)
722  {
723  LhsPacket A0, A1;
724  RhsPacket B_0;
725  RhsPacket T0;
726 
727  traits.loadLhs(&blA[0*LhsProgress], A0);
728  traits.loadLhs(&blA[1*LhsProgress], A1);
729  traits.loadRhs(&blB[0*RhsProgress], B_0);
730  traits.madd(A0,B_0,C0,T0);
731  traits.madd(A1,B_0,C4,B_0);
732  traits.loadRhs(&blB[1*RhsProgress], B_0);
733  traits.madd(A0,B_0,C1,T0);
734  traits.madd(A1,B_0,C5,B_0);
735  }
736  else
737  {
738  LhsPacket A0, A1;
739  RhsPacket B_0, B1, B2, B3;
740  RhsPacket T0;
741 
742  traits.loadLhs(&blA[0*LhsProgress], A0);
743  traits.loadLhs(&blA[1*LhsProgress], A1);
744  traits.loadRhs(&blB[0*RhsProgress], B_0);
745  traits.loadRhs(&blB[1*RhsProgress], B1);
746 
747  traits.madd(A0,B_0,C0,T0);
748  traits.loadRhs(&blB[2*RhsProgress], B2);
749  traits.madd(A1,B_0,C4,B_0);
750  traits.loadRhs(&blB[3*RhsProgress], B3);
751  traits.madd(A0,B1,C1,T0);
752  traits.madd(A1,B1,C5,B1);
753  traits.madd(A0,B2,C2,T0);
754  traits.madd(A1,B2,C6,B2);
755  traits.madd(A0,B3,C3,T0);
756  traits.madd(A1,B3,C7,B3);
757  }
758 
759  blB += nr*RhsProgress;
760  blA += mr;
761  }
762 
763  if(nr==4)
764  {
765  ResPacket R0, R1, R2, R3, R4, R5, R6;
766  ResPacket alphav = pset1<ResPacket>(alpha);
767 
768  R0 = ploadu<ResPacket>(r0);
769  R1 = ploadu<ResPacket>(r1);
770  R2 = ploadu<ResPacket>(r2);
771  R3 = ploadu<ResPacket>(r3);
772  R4 = ploadu<ResPacket>(r0 + ResPacketSize);
773  R5 = ploadu<ResPacket>(r1 + ResPacketSize);
774  R6 = ploadu<ResPacket>(r2 + ResPacketSize);
775  traits.acc(C0, alphav, R0);
776  pstoreu(r0, R0);
777  R0 = ploadu<ResPacket>(r3 + ResPacketSize);
778 
779  traits.acc(C1, alphav, R1);
780  traits.acc(C2, alphav, R2);
781  traits.acc(C3, alphav, R3);
782  traits.acc(C4, alphav, R4);
783  traits.acc(C5, alphav, R5);
784  traits.acc(C6, alphav, R6);
785  traits.acc(C7, alphav, R0);
786 
787  pstoreu(r1, R1);
788  pstoreu(r2, R2);
789  pstoreu(r3, R3);
790  pstoreu(r0 + ResPacketSize, R4);
791  pstoreu(r1 + ResPacketSize, R5);
792  pstoreu(r2 + ResPacketSize, R6);
793  pstoreu(r3 + ResPacketSize, R0);
794  }
795  else
796  {
797  ResPacket R0, R1, R4;
798  ResPacket alphav = pset1<ResPacket>(alpha);
799 
800  R0 = ploadu<ResPacket>(r0);
801  R1 = ploadu<ResPacket>(r1);
802  R4 = ploadu<ResPacket>(r0 + ResPacketSize);
803  traits.acc(C0, alphav, R0);
804  pstoreu(r0, R0);
805  R0 = ploadu<ResPacket>(r1 + ResPacketSize);
806  traits.acc(C1, alphav, R1);
807  traits.acc(C4, alphav, R4);
808  traits.acc(C5, alphav, R0);
809  pstoreu(r1, R1);
810  pstoreu(r0 + ResPacketSize, R4);
811  pstoreu(r1 + ResPacketSize, R0);
812  }
813 
814  }
815 
816  if(rows-peeled_mc>=LhsProgress)
817  {
818  Index i = peeled_mc;
819  const LhsScalar* blA = &blockA[i*strideA+offsetA*LhsProgress];
820  prefetch(&blA[0]);
821 
822  // gets res block as register
823  AccPacket C0, C1, C2, C3;
824  traits.initAcc(C0);
825  traits.initAcc(C1);
826  if(nr==4) traits.initAcc(C2);
827  if(nr==4) traits.initAcc(C3);
828 
829  // performs "inner" product
830  const RhsScalar* blB = unpackedB;
831  for(Index k=0; k<peeled_kc; k+=4)
832  {
833  if(nr==2)
834  {
835  LhsPacket A0;
836  RhsPacket B_0, B1;
837 
838  traits.loadLhs(&blA[0*LhsProgress], A0);
839  traits.loadRhs(&blB[0*RhsProgress], B_0);
840  traits.loadRhs(&blB[1*RhsProgress], B1);
841  traits.madd(A0,B_0,C0,B_0);
842  traits.loadRhs(&blB[2*RhsProgress], B_0);
843  traits.madd(A0,B1,C1,B1);
844  traits.loadLhs(&blA[1*LhsProgress], A0);
845  traits.loadRhs(&blB[3*RhsProgress], B1);
846  traits.madd(A0,B_0,C0,B_0);
847  traits.loadRhs(&blB[4*RhsProgress], B_0);
848  traits.madd(A0,B1,C1,B1);
849  traits.loadLhs(&blA[2*LhsProgress], A0);
850  traits.loadRhs(&blB[5*RhsProgress], B1);
851  traits.madd(A0,B_0,C0,B_0);
852  traits.loadRhs(&blB[6*RhsProgress], B_0);
853  traits.madd(A0,B1,C1,B1);
854  traits.loadLhs(&blA[3*LhsProgress], A0);
855  traits.loadRhs(&blB[7*RhsProgress], B1);
856  traits.madd(A0,B_0,C0,B_0);
857  traits.madd(A0,B1,C1,B1);
858  }
859  else
860  {
861  LhsPacket A0;
862  RhsPacket B_0, B1, B2, B3;
863 
864  traits.loadLhs(&blA[0*LhsProgress], A0);
865  traits.loadRhs(&blB[0*RhsProgress], B_0);
866  traits.loadRhs(&blB[1*RhsProgress], B1);
867 
868  traits.madd(A0,B_0,C0,B_0);
869  traits.loadRhs(&blB[2*RhsProgress], B2);
870  traits.loadRhs(&blB[3*RhsProgress], B3);
871  traits.loadRhs(&blB[4*RhsProgress], B_0);
872  traits.madd(A0,B1,C1,B1);
873  traits.loadRhs(&blB[5*RhsProgress], B1);
874  traits.madd(A0,B2,C2,B2);
875  traits.loadRhs(&blB[6*RhsProgress], B2);
876  traits.madd(A0,B3,C3,B3);
877  traits.loadLhs(&blA[1*LhsProgress], A0);
878  traits.loadRhs(&blB[7*RhsProgress], B3);
879  traits.madd(A0,B_0,C0,B_0);
880  traits.loadRhs(&blB[8*RhsProgress], B_0);
881  traits.madd(A0,B1,C1,B1);
882  traits.loadRhs(&blB[9*RhsProgress], B1);
883  traits.madd(A0,B2,C2,B2);
884  traits.loadRhs(&blB[10*RhsProgress], B2);
885  traits.madd(A0,B3,C3,B3);
886  traits.loadLhs(&blA[2*LhsProgress], A0);
887  traits.loadRhs(&blB[11*RhsProgress], B3);
888 
889  traits.madd(A0,B_0,C0,B_0);
890  traits.loadRhs(&blB[12*RhsProgress], B_0);
891  traits.madd(A0,B1,C1,B1);
892  traits.loadRhs(&blB[13*RhsProgress], B1);
893  traits.madd(A0,B2,C2,B2);
894  traits.loadRhs(&blB[14*RhsProgress], B2);
895  traits.madd(A0,B3,C3,B3);
896 
897  traits.loadLhs(&blA[3*LhsProgress], A0);
898  traits.loadRhs(&blB[15*RhsProgress], B3);
899  traits.madd(A0,B_0,C0,B_0);
900  traits.madd(A0,B1,C1,B1);
901  traits.madd(A0,B2,C2,B2);
902  traits.madd(A0,B3,C3,B3);
903  }
904 
905  blB += nr*4*RhsProgress;
906  blA += 4*LhsProgress;
907  }
908  // process remaining peeled loop
909  for(Index k=peeled_kc; k<depth; k++)
910  {
911  if(nr==2)
912  {
913  LhsPacket A0;
914  RhsPacket B_0, B1;
915 
916  traits.loadLhs(&blA[0*LhsProgress], A0);
917  traits.loadRhs(&blB[0*RhsProgress], B_0);
918  traits.loadRhs(&blB[1*RhsProgress], B1);
919  traits.madd(A0,B_0,C0,B_0);
920  traits.madd(A0,B1,C1,B1);
921  }
922  else
923  {
924  LhsPacket A0;
925  RhsPacket B_0, B1, B2, B3;
926 
927  traits.loadLhs(&blA[0*LhsProgress], A0);
928  traits.loadRhs(&blB[0*RhsProgress], B_0);
929  traits.loadRhs(&blB[1*RhsProgress], B1);
930  traits.loadRhs(&blB[2*RhsProgress], B2);
931  traits.loadRhs(&blB[3*RhsProgress], B3);
932 
933  traits.madd(A0,B_0,C0,B_0);
934  traits.madd(A0,B1,C1,B1);
935  traits.madd(A0,B2,C2,B2);
936  traits.madd(A0,B3,C3,B3);
937  }
938 
939  blB += nr*RhsProgress;
940  blA += LhsProgress;
941  }
942 
943  ResPacket R0, R1, R2, R3;
944  ResPacket alphav = pset1<ResPacket>(alpha);
945 
946  ResScalar* r0 = &res[(j2+0)*resStride + i];
947  ResScalar* r1 = r0 + resStride;
948  ResScalar* r2 = r1 + resStride;
949  ResScalar* r3 = r2 + resStride;
950 
951  R0 = ploadu<ResPacket>(r0);
952  R1 = ploadu<ResPacket>(r1);
953  if(nr==4) R2 = ploadu<ResPacket>(r2);
954  if(nr==4) R3 = ploadu<ResPacket>(r3);
955 
956  traits.acc(C0, alphav, R0);
957  traits.acc(C1, alphav, R1);
958  if(nr==4) traits.acc(C2, alphav, R2);
959  if(nr==4) traits.acc(C3, alphav, R3);
960 
961  pstoreu(r0, R0);
962  pstoreu(r1, R1);
963  if(nr==4) pstoreu(r2, R2);
964  if(nr==4) pstoreu(r3, R3);
965  }
966  for(Index i=peeled_mc2; i<rows; i++)
967  {
968  const LhsScalar* blA = &blockA[i*strideA+offsetA];
969  prefetch(&blA[0]);
970 
971  // gets a 1 x nr res block as registers
972  ResScalar C0(0), C1(0), C2(0), C3(0);
973  // TODO directly use blockB ???
974  const RhsScalar* blB = &blockB[j2*strideB+offsetB*nr];
975  for(Index k=0; k<depth; k++)
976  {
977  if(nr==2)
978  {
979  LhsScalar A0;
980  RhsScalar B_0, B1;
981 
982  A0 = blA[k];
983  B_0 = blB[0];
984  B1 = blB[1];
985  MADD(cj,A0,B_0,C0,B_0);
986  MADD(cj,A0,B1,C1,B1);
987  }
988  else
989  {
990  LhsScalar A0;
991  RhsScalar B_0, B1, B2, B3;
992 
993  A0 = blA[k];
994  B_0 = blB[0];
995  B1 = blB[1];
996  B2 = blB[2];
997  B3 = blB[3];
998 
999  MADD(cj,A0,B_0,C0,B_0);
1000  MADD(cj,A0,B1,C1,B1);
1001  MADD(cj,A0,B2,C2,B2);
1002  MADD(cj,A0,B3,C3,B3);
1003  }
1004 
1005  blB += nr;
1006  }
1007  res[(j2+0)*resStride + i] += alpha*C0;
1008  res[(j2+1)*resStride + i] += alpha*C1;
1009  if(nr==4) res[(j2+2)*resStride + i] += alpha*C2;
1010  if(nr==4) res[(j2+3)*resStride + i] += alpha*C3;
1011  }
1012  }
1013  // process remaining rhs/res columns one at a time
1014  // => do the same but with nr==1
1015  for(Index j2=packet_cols; j2<cols; j2++)
1016  {
1017  // unpack B
1018  traits.unpackRhs(depth, &blockB[j2*strideB+offsetB], unpackedB);
1019 
1020  for(Index i=0; i<peeled_mc; i+=mr)
1021  {
1022  const LhsScalar* blA = &blockA[i*strideA+offsetA*mr];
1023  prefetch(&blA[0]);
1024 
1025  // TODO move the res loads to the stores
1026 
1027  // get res block as registers
1028  AccPacket C0, C4;
1029  traits.initAcc(C0);
1030  traits.initAcc(C4);
1031 
1032  const RhsScalar* blB = unpackedB;
1033  for(Index k=0; k<depth; k++)
1034  {
1035  LhsPacket A0, A1;
1036  RhsPacket B_0;
1037  RhsPacket T0;
1038 
1039  traits.loadLhs(&blA[0*LhsProgress], A0);
1040  traits.loadLhs(&blA[1*LhsProgress], A1);
1041  traits.loadRhs(&blB[0*RhsProgress], B_0);
1042  traits.madd(A0,B_0,C0,T0);
1043  traits.madd(A1,B_0,C4,B_0);
1044 
1045  blB += RhsProgress;
1046  blA += 2*LhsProgress;
1047  }
1048  ResPacket R0, R4;
1049  ResPacket alphav = pset1<ResPacket>(alpha);
1050 
1051  ResScalar* r0 = &res[(j2+0)*resStride + i];
1052 
1053  R0 = ploadu<ResPacket>(r0);
1054  R4 = ploadu<ResPacket>(r0+ResPacketSize);
1055 
1056  traits.acc(C0, alphav, R0);
1057  traits.acc(C4, alphav, R4);
1058 
1059  pstoreu(r0, R0);
1060  pstoreu(r0+ResPacketSize, R4);
1061  }
1062  if(rows-peeled_mc>=LhsProgress)
1063  {
1064  Index i = peeled_mc;
1065  const LhsScalar* blA = &blockA[i*strideA+offsetA*LhsProgress];
1066  prefetch(&blA[0]);
1067 
1068  AccPacket C0;
1069  traits.initAcc(C0);
1070 
1071  const RhsScalar* blB = unpackedB;
1072  for(Index k=0; k<depth; k++)
1073  {
1074  LhsPacket A0;
1075  RhsPacket B_0;
1076  traits.loadLhs(blA, A0);
1077  traits.loadRhs(blB, B_0);
1078  traits.madd(A0, B_0, C0, B_0);
1079  blB += RhsProgress;
1080  blA += LhsProgress;
1081  }
1082 
1083  ResPacket alphav = pset1<ResPacket>(alpha);
1084  ResPacket R0 = ploadu<ResPacket>(&res[(j2+0)*resStride + i]);
1085  traits.acc(C0, alphav, R0);
1086  pstoreu(&res[(j2+0)*resStride + i], R0);
1087  }
1088  for(Index i=peeled_mc2; i<rows; i++)
1089  {
1090  const LhsScalar* blA = &blockA[i*strideA+offsetA];
1091  prefetch(&blA[0]);
1092 
1093  // gets a 1 x 1 res block as registers
1094  ResScalar C0(0);
1095  // FIXME directly use blockB ??
1096  const RhsScalar* blB = &blockB[j2*strideB+offsetB];
1097  for(Index k=0; k<depth; k++)
1098  {
1099  LhsScalar A0 = blA[k];
1100  RhsScalar B_0 = blB[k];
1101  MADD(cj, A0, B_0, C0, B_0);
1102  }
1103  res[(j2+0)*resStride + i] += alpha*C0;
1104  }
1105  }
1106  }
1107 };
1108 
1109 #undef CJMADD
1110 
1111 // pack a block of the lhs
1112 // The traversal is as follow (mr==4):
1113 // 0 4 8 12 ...
1114 // 1 5 9 13 ...
1115 // 2 6 10 14 ...
1116 // 3 7 11 15 ...
1117 //
1118 // 16 20 24 28 ...
1119 // 17 21 25 29 ...
1120 // 18 22 26 30 ...
1121 // 19 23 27 31 ...
1122 //
1123 // 32 33 34 35 ...
1124 // 36 36 38 39 ...
1125 template<typename Scalar, typename Index, int Pack1, int Pack2, int StorageOrder, bool Conjugate, bool PanelMode>
1126 struct gemm_pack_lhs
1127 {
1128  EIGEN_DONT_INLINE void operator()(Scalar* blockA, const Scalar* EIGEN_RESTRICT _lhs, Index lhsStride, Index depth, Index rows,
1129  Index stride=0, Index offset=0)
1130  {
1131  typedef typename packet_traits<Scalar>::type Packet;
1132  enum { PacketSize = packet_traits<Scalar>::size };
1133 
1134  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS");
1135  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1136  eigen_assert( (StorageOrder==RowMajor) || ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) );
1137  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1138  const_blas_data_mapper<Scalar, Index, StorageOrder> lhs(_lhs,lhsStride);
1139  Index count = 0;
1140  Index peeled_mc = (rows/Pack1)*Pack1;
1141  for(Index i=0; i<peeled_mc; i+=Pack1)
1142  {
1143  if(PanelMode) count += Pack1 * offset;
1144 
1145  if(StorageOrder==ColMajor)
1146  {
1147  for(Index k=0; k<depth; k++)
1148  {
1149  Packet A, B, C, D;
1150  if(Pack1>=1*PacketSize) A = ploadu<Packet>(&lhs(i+0*PacketSize, k));
1151  if(Pack1>=2*PacketSize) B = ploadu<Packet>(&lhs(i+1*PacketSize, k));
1152  if(Pack1>=3*PacketSize) C = ploadu<Packet>(&lhs(i+2*PacketSize, k));
1153  if(Pack1>=4*PacketSize) D = ploadu<Packet>(&lhs(i+3*PacketSize, k));
1154  if(Pack1>=1*PacketSize) { pstore(blockA+count, cj.pconj(A)); count+=PacketSize; }
1155  if(Pack1>=2*PacketSize) { pstore(blockA+count, cj.pconj(B)); count+=PacketSize; }
1156  if(Pack1>=3*PacketSize) { pstore(blockA+count, cj.pconj(C)); count+=PacketSize; }
1157  if(Pack1>=4*PacketSize) { pstore(blockA+count, cj.pconj(D)); count+=PacketSize; }
1158  }
1159  }
1160  else
1161  {
1162  for(Index k=0; k<depth; k++)
1163  {
1164  // TODO add a vectorized transpose here
1165  Index w=0;
1166  for(; w<Pack1-3; w+=4)
1167  {
1168  Scalar a(cj(lhs(i+w+0, k))),
1169  b(cj(lhs(i+w+1, k))),
1170  c(cj(lhs(i+w+2, k))),
1171  d(cj(lhs(i+w+3, k)));
1172  blockA[count++] = a;
1173  blockA[count++] = b;
1174  blockA[count++] = c;
1175  blockA[count++] = d;
1176  }
1177  if(Pack1%4)
1178  for(;w<Pack1;++w)
1179  blockA[count++] = cj(lhs(i+w, k));
1180  }
1181  }
1182  if(PanelMode) count += Pack1 * (stride-offset-depth);
1183  }
1184  if(rows-peeled_mc>=Pack2)
1185  {
1186  if(PanelMode) count += Pack2*offset;
1187  for(Index k=0; k<depth; k++)
1188  for(Index w=0; w<Pack2; w++)
1189  blockA[count++] = cj(lhs(peeled_mc+w, k));
1190  if(PanelMode) count += Pack2 * (stride-offset-depth);
1191  peeled_mc += Pack2;
1192  }
1193  for(Index i=peeled_mc; i<rows; i++)
1194  {
1195  if(PanelMode) count += offset;
1196  for(Index k=0; k<depth; k++)
1197  blockA[count++] = cj(lhs(i, k));
1198  if(PanelMode) count += (stride-offset-depth);
1199  }
1200  }
1201 };
1202 
1203 // copy a complete panel of the rhs
1204 // this version is optimized for column major matrices
1205 // The traversal order is as follow: (nr==4):
1206 // 0 1 2 3 12 13 14 15 24 27
1207 // 4 5 6 7 16 17 18 19 25 28
1208 // 8 9 10 11 20 21 22 23 26 29
1209 // . . . . . . . . . .
1210 template<typename Scalar, typename Index, int nr, bool Conjugate, bool PanelMode>
1211 struct gemm_pack_rhs<Scalar, Index, nr, ColMajor, Conjugate, PanelMode>
1212 {
1213  typedef typename packet_traits<Scalar>::type Packet;
1214  enum { PacketSize = packet_traits<Scalar>::size };
1215  EIGEN_DONT_INLINE void operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols,
1216  Index stride=0, Index offset=0)
1217  {
1218  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS COLMAJOR");
1219  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1220  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1221  Index packet_cols = (cols/nr) * nr;
1222  Index count = 0;
1223  for(Index j2=0; j2<packet_cols; j2+=nr)
1224  {
1225  // skip what we have before
1226  if(PanelMode) count += nr * offset;
1227  const Scalar* b0 = &rhs[(j2+0)*rhsStride];
1228  const Scalar* b1 = &rhs[(j2+1)*rhsStride];
1229  const Scalar* b2 = &rhs[(j2+2)*rhsStride];
1230  const Scalar* b3 = &rhs[(j2+3)*rhsStride];
1231  for(Index k=0; k<depth; k++)
1232  {
1233  blockB[count+0] = cj(b0[k]);
1234  blockB[count+1] = cj(b1[k]);
1235  if(nr==4) blockB[count+2] = cj(b2[k]);
1236  if(nr==4) blockB[count+3] = cj(b3[k]);
1237  count += nr;
1238  }
1239  // skip what we have after
1240  if(PanelMode) count += nr * (stride-offset-depth);
1241  }
1242 
1243  // copy the remaining columns one at a time (nr==1)
1244  for(Index j2=packet_cols; j2<cols; ++j2)
1245  {
1246  if(PanelMode) count += offset;
1247  const Scalar* b0 = &rhs[(j2+0)*rhsStride];
1248  for(Index k=0; k<depth; k++)
1249  {
1250  blockB[count] = cj(b0[k]);
1251  count += 1;
1252  }
1253  if(PanelMode) count += (stride-offset-depth);
1254  }
1255  }
1256 };
1257 
1258 // this version is optimized for row major matrices
1259 template<typename Scalar, typename Index, int nr, bool Conjugate, bool PanelMode>
1260 struct gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, PanelMode>
1261 {
1262  enum { PacketSize = packet_traits<Scalar>::size };
1263  EIGEN_DONT_INLINE void operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols,
1264  Index stride=0, Index offset=0)
1265  {
1266  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS ROWMAJOR");
1267  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1268  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1269  Index packet_cols = (cols/nr) * nr;
1270  Index count = 0;
1271  for(Index j2=0; j2<packet_cols; j2+=nr)
1272  {
1273  // skip what we have before
1274  if(PanelMode) count += nr * offset;
1275  for(Index k=0; k<depth; k++)
1276  {
1277  const Scalar* b0 = &rhs[k*rhsStride + j2];
1278  blockB[count+0] = cj(b0[0]);
1279  blockB[count+1] = cj(b0[1]);
1280  if(nr==4) blockB[count+2] = cj(b0[2]);
1281  if(nr==4) blockB[count+3] = cj(b0[3]);
1282  count += nr;
1283  }
1284  // skip what we have after
1285  if(PanelMode) count += nr * (stride-offset-depth);
1286  }
1287  // copy the remaining columns one at a time (nr==1)
1288  for(Index j2=packet_cols; j2<cols; ++j2)
1289  {
1290  if(PanelMode) count += offset;
1291  const Scalar* b0 = &rhs[j2];
1292  for(Index k=0; k<depth; k++)
1293  {
1294  blockB[count] = cj(b0[k*rhsStride]);
1295  count += 1;
1296  }
1297  if(PanelMode) count += stride-offset-depth;
1298  }
1299  }
1300 };
1301 
1302 } // end namespace internal
1303 
1306 inline std::ptrdiff_t l1CacheSize()
1307 {
1308  std::ptrdiff_t l1, l2;
1310  return l1;
1311 }
1312 
1315 inline std::ptrdiff_t l2CacheSize()
1316 {
1317  std::ptrdiff_t l1, l2;
1319  return l2;
1320 }
1321 
1327 inline void setCpuCacheSizes(std::ptrdiff_t l1, std::ptrdiff_t l2)
1328 {
1330 }
1331 
1332 } // end namespace Eigen
1333 
1334 #endif // EIGEN_GENERAL_BLOCK_PANEL_H