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 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_GENERAL_BLOCK_PANEL_H
11 #define EIGEN_GENERAL_BLOCK_PANEL_H
12 
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
18 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs=false, bool _ConjRhs=false>
19 class gebp_traits;
20 
21 
23 inline std::ptrdiff_t manage_caching_sizes_helper(std::ptrdiff_t a, std::ptrdiff_t b)
24 {
25  return a<=0 ? b : a;
26 }
27 
28 #if EIGEN_ARCH_i386_OR_x86_64
29 const std::ptrdiff_t defaultL1CacheSize = 32*1024;
30 const std::ptrdiff_t defaultL2CacheSize = 256*1024;
31 const std::ptrdiff_t defaultL3CacheSize = 2*1024*1024;
32 #else
33 const std::ptrdiff_t defaultL1CacheSize = 16*1024;
34 const std::ptrdiff_t defaultL2CacheSize = 512*1024;
35 const std::ptrdiff_t defaultL3CacheSize = 512*1024;
36 #endif
37 
39 struct CacheSizes {
40  CacheSizes(): m_l1(-1),m_l2(-1),m_l3(-1) {
41  int l1CacheSize, l2CacheSize, l3CacheSize;
42  queryCacheSizes(l1CacheSize, l2CacheSize, l3CacheSize);
43  m_l1 = manage_caching_sizes_helper(l1CacheSize, defaultL1CacheSize);
44  m_l2 = manage_caching_sizes_helper(l2CacheSize, defaultL2CacheSize);
45  m_l3 = manage_caching_sizes_helper(l3CacheSize, defaultL3CacheSize);
46  }
47 
48  std::ptrdiff_t m_l1;
49  std::ptrdiff_t m_l2;
50  std::ptrdiff_t m_l3;
51 };
52 
53 
55 inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1, std::ptrdiff_t* l2, std::ptrdiff_t* l3)
56 {
57  static CacheSizes m_cacheSizes;
58 
59  if(action==SetAction)
60  {
61  // set the cpu cache size and cache all block sizes from a global cache size in byte
62  eigen_internal_assert(l1!=0 && l2!=0);
63  m_cacheSizes.m_l1 = *l1;
64  m_cacheSizes.m_l2 = *l2;
65  m_cacheSizes.m_l3 = *l3;
66  }
67  else if(action==GetAction)
68  {
69  eigen_internal_assert(l1!=0 && l2!=0);
70  *l1 = m_cacheSizes.m_l1;
71  *l2 = m_cacheSizes.m_l2;
72  *l3 = m_cacheSizes.m_l3;
73  }
74  else
75  {
76  eigen_internal_assert(false);
77  }
78 }
79 
80 /* Helper for computeProductBlockingSizes.
81  *
82  * Given a m x k times k x n matrix product of scalar types \c LhsScalar and \c RhsScalar,
83  * this function computes the blocking size parameters along the respective dimensions
84  * for matrix products and related algorithms. The blocking sizes depends on various
85  * parameters:
86  * - the L1 and L2 cache sizes,
87  * - the register level blocking sizes defined by gebp_traits,
88  * - the number of scalars that fit into a packet (when vectorization is enabled).
89  *
90  * \sa setCpuCacheSizes */
91 
92 template<typename LhsScalar, typename RhsScalar, int KcFactor, typename Index>
93 void evaluateProductBlockingSizesHeuristic(Index& k, Index& m, Index& n, Index num_threads = 1)
94 {
95  typedef gebp_traits<LhsScalar,RhsScalar> Traits;
96 
97  // Explanations:
98  // Let's recall that the product algorithms form mc x kc vertical panels A' on the lhs and
99  // kc x nc blocks B' on the rhs. B' has to fit into L2/L3 cache. Moreover, A' is processed
100  // per mr x kc horizontal small panels where mr is the blocking size along the m dimension
101  // at the register level. This small horizontal panel has to stay within L1 cache.
102  std::ptrdiff_t l1, l2, l3;
103  manage_caching_sizes(GetAction, &l1, &l2, &l3);
104 
105  if (num_threads > 1) {
106  typedef typename Traits::ResScalar ResScalar;
107  enum {
108  kdiv = KcFactor * (Traits::mr * sizeof(LhsScalar) + Traits::nr * sizeof(RhsScalar)),
109  ksub = Traits::mr * Traits::nr * sizeof(ResScalar),
110  kr = 8,
111  mr = Traits::mr,
112  nr = Traits::nr
113  };
114  // Increasing k gives us more time to prefetch the content of the "C"
115  // registers. However once the latency is hidden there is no point in
116  // increasing the value of k, so we'll cap it at 320 (value determined
117  // experimentally).
118  const Index k_cache = (numext::mini<Index>)((l1-ksub)/kdiv, 320);
119  if (k_cache < k) {
120  k = k_cache - (k_cache % kr);
121  eigen_internal_assert(k > 0);
122  }
123 
124  const Index n_cache = (l2-l1) / (nr * sizeof(RhsScalar) * k);
125  const Index n_per_thread = numext::div_ceil(n, num_threads);
126  if (n_cache <= n_per_thread) {
127  // Don't exceed the capacity of the l2 cache.
128  eigen_internal_assert(n_cache >= static_cast<Index>(nr));
129  n = n_cache - (n_cache % nr);
130  eigen_internal_assert(n > 0);
131  } else {
132  n = (numext::mini<Index>)(n, (n_per_thread + nr - 1) - ((n_per_thread + nr - 1) % nr));
133  }
134 
135  if (l3 > l2) {
136  // l3 is shared between all cores, so we'll give each thread its own chunk of l3.
137  const Index m_cache = (l3-l2) / (sizeof(LhsScalar) * k * num_threads);
138  const Index m_per_thread = numext::div_ceil(m, num_threads);
139  if(m_cache < m_per_thread && m_cache >= static_cast<Index>(mr)) {
140  m = m_cache - (m_cache % mr);
141  eigen_internal_assert(m > 0);
142  } else {
143  m = (numext::mini<Index>)(m, (m_per_thread + mr - 1) - ((m_per_thread + mr - 1) % mr));
144  }
145  }
146  }
147  else {
148  // In unit tests we do not want to use extra large matrices,
149  // so we reduce the cache size to check the blocking strategy is not flawed
150 #ifdef EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS
151  l1 = 9*1024;
152  l2 = 32*1024;
153  l3 = 512*1024;
154 #endif
155 
156  // Early return for small problems because the computation below are time consuming for small problems.
157  // Perhaps it would make more sense to consider k*n*m??
158  // Note that for very tiny problem, this function should be bypassed anyway
159  // because we use the coefficient-based implementation for them.
160  if((numext::maxi)(k,(numext::maxi)(m,n))<48)
161  return;
162 
163  typedef typename Traits::ResScalar ResScalar;
164  enum {
165  k_peeling = 8,
166  k_div = KcFactor * (Traits::mr * sizeof(LhsScalar) + Traits::nr * sizeof(RhsScalar)),
167  k_sub = Traits::mr * Traits::nr * sizeof(ResScalar)
168  };
169 
170  // ---- 1st level of blocking on L1, yields kc ----
171 
172  // Blocking on the third dimension (i.e., k) is chosen so that an horizontal panel
173  // of size mr x kc of the lhs plus a vertical panel of kc x nr of the rhs both fits within L1 cache.
174  // We also include a register-level block of the result (mx x nr).
175  // (In an ideal world only the lhs panel would stay in L1)
176  // Moreover, kc has to be a multiple of 8 to be compatible with loop peeling, leading to a maximum blocking size of:
177  const Index max_kc = numext::maxi<Index>(((l1-k_sub)/k_div) & (~(k_peeling-1)),1);
178  const Index old_k = k;
179  if(k>max_kc)
180  {
181  // We are really blocking on the third dimension:
182  // -> reduce blocking size to make sure the last block is as large as possible
183  // while keeping the same number of sweeps over the result.
184  k = (k%max_kc)==0 ? max_kc
185  : max_kc - k_peeling * ((max_kc-1-(k%max_kc))/(k_peeling*(k/max_kc+1)));
186 
187  eigen_internal_assert(((old_k/k) == (old_k/max_kc)) && "the number of sweeps has to remain the same");
188  }
189 
190  // ---- 2nd level of blocking on max(L2,L3), yields nc ----
191 
192  // TODO find a reliable way to get the actual amount of cache per core to use for 2nd level blocking, that is:
193  // actual_l2 = max(l2, l3/nb_core_sharing_l3)
194  // The number below is quite conservative: it is better to underestimate the cache size rather than overestimating it)
195  // For instance, it corresponds to 6MB of L3 shared among 4 cores.
196  #ifdef EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS
197  const Index actual_l2 = l3;
198  #else
199  const Index actual_l2 = 1572864; // == 1.5 MB
200  #endif
201 
202  // Here, nc is chosen such that a block of kc x nc of the rhs fit within half of L2.
203  // The second half is implicitly reserved to access the result and lhs coefficients.
204  // When k<max_kc, then nc can arbitrarily growth. In practice, it seems to be fruitful
205  // to limit this growth: we bound nc to growth by a factor x1.5.
206  // However, if the entire lhs block fit within L1, then we are not going to block on the rows at all,
207  // and it becomes fruitful to keep the packed rhs blocks in L1 if there is enough remaining space.
208  Index max_nc;
209  const Index lhs_bytes = m * k * sizeof(LhsScalar);
210  const Index remaining_l1 = l1- k_sub - lhs_bytes;
211  if(remaining_l1 >= Index(Traits::nr*sizeof(RhsScalar))*k)
212  {
213  // L1 blocking
214  max_nc = remaining_l1 / (k*sizeof(RhsScalar));
215  }
216  else
217  {
218  // L2 blocking
219  max_nc = (3*actual_l2)/(2*2*max_kc*sizeof(RhsScalar));
220  }
221  // WARNING Below, we assume that Traits::nr is a power of two.
222  Index nc = numext::mini<Index>(actual_l2/(2*k*sizeof(RhsScalar)), max_nc) & (~(Traits::nr-1));
223  if(n>nc)
224  {
225  // We are really blocking over the columns:
226  // -> reduce blocking size to make sure the last block is as large as possible
227  // while keeping the same number of sweeps over the packed lhs.
228  // Here we allow one more sweep if this gives us a perfect match, thus the commented "-1"
229  n = (n%nc)==0 ? nc
230  : (nc - Traits::nr * ((nc/*-1*/-(n%nc))/(Traits::nr*(n/nc+1))));
231  }
232  else if(old_k==k)
233  {
234  // So far, no blocking at all, i.e., kc==k, and nc==n.
235  // In this case, let's perform a blocking over the rows such that the packed lhs data is kept in cache L1/L2
236  // TODO: part of this blocking strategy is now implemented within the kernel itself, so the L1-based heuristic here should be obsolete.
237  Index problem_size = k*n*sizeof(LhsScalar);
238  Index actual_lm = actual_l2;
239  Index max_mc = m;
240  if(problem_size<=1024)
241  {
242  // problem is small enough to keep in L1
243  // Let's choose m such that lhs's block fit in 1/3 of L1
244  actual_lm = l1;
245  }
246  else if(l3!=0 && problem_size<=32768)
247  {
248  // we have both L2 and L3, and problem is small enough to be kept in L2
249  // Let's choose m such that lhs's block fit in 1/3 of L2
250  actual_lm = l2;
251  max_mc = (numext::mini<Index>)(576,max_mc);
252  }
253  Index mc = (numext::mini<Index>)(actual_lm/(3*k*sizeof(LhsScalar)), max_mc);
254  if (mc > Traits::mr) mc -= mc % Traits::mr;
255  else if (mc==0) return;
256  m = (m%mc)==0 ? mc
257  : (mc - Traits::mr * ((mc/*-1*/-(m%mc))/(Traits::mr*(m/mc+1))));
258  }
259  }
260 }
261 
262 template <typename Index>
263 inline bool useSpecificBlockingSizes(Index& k, Index& m, Index& n)
264 {
265 #ifdef EIGEN_TEST_SPECIFIC_BLOCKING_SIZES
266  if (EIGEN_TEST_SPECIFIC_BLOCKING_SIZES) {
267  k = numext::mini<Index>(k, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_K);
268  m = numext::mini<Index>(m, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_M);
269  n = numext::mini<Index>(n, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_N);
270  return true;
271  }
272 #else
273  EIGEN_UNUSED_VARIABLE(k)
274  EIGEN_UNUSED_VARIABLE(m)
275  EIGEN_UNUSED_VARIABLE(n)
276 #endif
277  return false;
278 }
279 
296 template<typename LhsScalar, typename RhsScalar, int KcFactor, typename Index>
297 void computeProductBlockingSizes(Index& k, Index& m, Index& n, Index num_threads = 1)
298 {
299  if (!useSpecificBlockingSizes(k, m, n)) {
300  evaluateProductBlockingSizesHeuristic<LhsScalar, RhsScalar, KcFactor, Index>(k, m, n, num_threads);
301  }
302 }
303 
304 template<typename LhsScalar, typename RhsScalar, typename Index>
305 inline void computeProductBlockingSizes(Index& k, Index& m, Index& n, Index num_threads = 1)
306 {
307  computeProductBlockingSizes<LhsScalar,RhsScalar,1,Index>(k, m, n, num_threads);
308 }
309 
310 #ifdef EIGEN_HAS_SINGLE_INSTRUCTION_CJMADD
311  #define CJMADD(CJ,A,B,C,T) C = CJ.pmadd(A,B,C);
312 #else
313 
314  // FIXME (a bit overkill maybe ?)
315 
316  template<typename CJ, typename A, typename B, typename C, typename T> struct gebp_madd_selector {
317  EIGEN_ALWAYS_INLINE static void run(const CJ& cj, A& a, B& b, C& c, T& /*t*/)
318  {
319  c = cj.pmadd(a,b,c);
320  }
321  };
322 
323  template<typename CJ, typename T> struct gebp_madd_selector<CJ,T,T,T,T> {
324  EIGEN_ALWAYS_INLINE static void run(const CJ& cj, T& a, T& b, T& c, T& t)
325  {
326  t = b; t = cj.pmul(a,t); c = padd(c,t);
327  }
328  };
329 
330  template<typename CJ, typename A, typename B, typename C, typename T>
331  EIGEN_STRONG_INLINE void gebp_madd(const CJ& cj, A& a, B& b, C& c, T& t)
332  {
333  gebp_madd_selector<CJ,A,B,C,T>::run(cj,a,b,c,t);
334  }
335 
336  #define CJMADD(CJ,A,B,C,T) gebp_madd(CJ,A,B,C,T);
337 // #define CJMADD(CJ,A,B,C,T) T = B; T = CJ.pmul(A,T); C = padd(C,T);
338 #endif
339 
340 /* Vectorization logic
341  * real*real: unpack rhs to constant packets, ...
342  *
343  * 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),
344  * storing each res packet into two packets (2x2),
345  * at the end combine them: swap the second and addsub them
346  * cf*cf : same but with 2x4 blocks
347  * cplx*real : unpack rhs to constant packets, ...
348  * real*cplx : load lhs as (a0,a0,a1,a1), and mul as usual
349  */
350 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs, bool _ConjRhs>
351 class gebp_traits
352 {
353 public:
354  typedef _LhsScalar LhsScalar;
355  typedef _RhsScalar RhsScalar;
356  typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
357 
358  enum {
359  ConjLhs = _ConjLhs,
360  ConjRhs = _ConjRhs,
361  Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable,
362  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
363  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
364  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
365 
366  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
367 
368  // register block size along the N direction must be 1 or 4
369  nr = 4,
370 
371  // register block size along the M direction (currently, this one cannot be modified)
372  default_mr = (EIGEN_PLAIN_ENUM_MIN(16,NumberOfRegisters)/2/nr)*LhsPacketSize,
373 #if defined(EIGEN_HAS_SINGLE_INSTRUCTION_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX)
374  // we assume 16 registers
375  // See bug 992, if the scalar type is not vectorizable but that EIGEN_HAS_SINGLE_INSTRUCTION_MADD is defined,
376  // then using 3*LhsPacketSize triggers non-implemented paths in syrk.
377  mr = Vectorizable ? 3*LhsPacketSize : default_mr,
378 #else
379  mr = default_mr,
380 #endif
381 
382  LhsProgress = LhsPacketSize,
383  RhsProgress = 1
384  };
385 
386  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
387  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
388  typedef typename packet_traits<ResScalar>::type _ResPacket;
389 
390  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
391  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
392  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
393 
394  typedef ResPacket AccPacket;
395 
396  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
397  {
398  p = pset1<ResPacket>(ResScalar(0));
399  }
400 
401  EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1, RhsPacket& b2, RhsPacket& b3)
402  {
403  pbroadcast4(b, b0, b1, b2, b3);
404  }
405 
406 // EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1)
407 // {
408 // pbroadcast2(b, b0, b1);
409 // }
410 
411  template<typename RhsPacketType>
412  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const
413  {
414  dest = pset1<RhsPacketType>(*b);
415  }
416 
417  EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const
418  {
419  dest = ploadquad<RhsPacket>(b);
420  }
421 
422  template<typename LhsPacketType>
423  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacketType& dest) const
424  {
425  dest = pload<LhsPacketType>(a);
426  }
427 
428  template<typename LhsPacketType>
429  EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const
430  {
431  dest = ploadu<LhsPacketType>(a);
432  }
433 
434  template<typename LhsPacketType, typename RhsPacketType, typename AccPacketType>
435  EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, AccPacketType& tmp) const
436  {
437  conj_helper<LhsPacketType,RhsPacketType,ConjLhs,ConjRhs> cj;
438  // It would be a lot cleaner to call pmadd all the time. Unfortunately if we
439  // let gcc allocate the register in which to store the result of the pmul
440  // (in the case where there is no FMA) gcc fails to figure out how to avoid
441  // spilling register.
442 #ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
443  EIGEN_UNUSED_VARIABLE(tmp);
444  c = cj.pmadd(a,b,c);
445 #else
446  tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
447 #endif
448  }
449 
450  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
451  {
452  r = pmadd(c,alpha,r);
453  }
454 
455  template<typename ResPacketHalf>
456  EIGEN_STRONG_INLINE void acc(const ResPacketHalf& c, const ResPacketHalf& alpha, ResPacketHalf& r) const
457  {
458  r = pmadd(c,alpha,r);
459  }
460 
461 };
462 
463 template<typename RealScalar, bool _ConjLhs>
464 class gebp_traits<std::complex<RealScalar>, RealScalar, _ConjLhs, false>
465 {
466 public:
467  typedef std::complex<RealScalar> LhsScalar;
468  typedef RealScalar RhsScalar;
469  typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
470 
471  enum {
472  ConjLhs = _ConjLhs,
473  ConjRhs = false,
474  Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable,
475  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
476  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
477  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
478 
479  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
480  nr = 4,
481 #if defined(EIGEN_HAS_SINGLE_INSTRUCTION_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX)
482  // we assume 16 registers
483  mr = 3*LhsPacketSize,
484 #else
485  mr = (EIGEN_PLAIN_ENUM_MIN(16,NumberOfRegisters)/2/nr)*LhsPacketSize,
486 #endif
487 
488  LhsProgress = LhsPacketSize,
489  RhsProgress = 1
490  };
491 
492  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
493  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
494  typedef typename packet_traits<ResScalar>::type _ResPacket;
495 
496  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
497  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
498  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
499 
500  typedef ResPacket AccPacket;
501 
502  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
503  {
504  p = pset1<ResPacket>(ResScalar(0));
505  }
506 
507  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
508  {
509  dest = pset1<RhsPacket>(*b);
510  }
511 
512  EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const
513  {
514  dest = pset1<RhsPacket>(*b);
515  }
516 
517  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
518  {
519  dest = pload<LhsPacket>(a);
520  }
521 
522  EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacket& dest) const
523  {
524  dest = ploadu<LhsPacket>(a);
525  }
526 
527  EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1, RhsPacket& b2, RhsPacket& b3)
528  {
529  pbroadcast4(b, b0, b1, b2, b3);
530  }
531 
532 // EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1)
533 // {
534 // pbroadcast2(b, b0, b1);
535 // }
536 
537  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp) const
538  {
539  madd_impl(a, b, c, tmp, typename conditional<Vectorizable,true_type,false_type>::type());
540  }
541 
542  EIGEN_STRONG_INLINE void madd_impl(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp, const true_type&) const
543  {
544 #ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
545  EIGEN_UNUSED_VARIABLE(tmp);
546  c.v = pmadd(a.v,b,c.v);
547 #else
548  tmp = b; tmp = pmul(a.v,tmp); c.v = padd(c.v,tmp);
549 #endif
550  }
551 
552  EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/, const false_type&) const
553  {
554  c += a * b;
555  }
556 
557  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
558  {
559  r = cj.pmadd(c,alpha,r);
560  }
561 
562 protected:
563  conj_helper<ResPacket,ResPacket,ConjLhs,false> cj;
564 };
565 
566 template<typename Packet>
567 struct DoublePacket
568 {
569  Packet first;
570  Packet second;
571 };
572 
573 template<typename Packet>
574 DoublePacket<Packet> padd(const DoublePacket<Packet> &a, const DoublePacket<Packet> &b)
575 {
576  DoublePacket<Packet> res;
577  res.first = padd(a.first, b.first);
578  res.second = padd(a.second,b.second);
579  return res;
580 }
581 
582 template<typename Packet>
583 const DoublePacket<Packet>& predux_downto4(const DoublePacket<Packet> &a)
584 {
585  return a;
586 }
587 
588 template<typename Packet> struct unpacket_traits<DoublePacket<Packet> > { typedef DoublePacket<Packet> half; };
589 // template<typename Packet>
590 // DoublePacket<Packet> pmadd(const DoublePacket<Packet> &a, const DoublePacket<Packet> &b)
591 // {
592 // DoublePacket<Packet> res;
593 // res.first = padd(a.first, b.first);
594 // res.second = padd(a.second,b.second);
595 // return res;
596 // }
597 
598 template<typename RealScalar, bool _ConjLhs, bool _ConjRhs>
599 class gebp_traits<std::complex<RealScalar>, std::complex<RealScalar>, _ConjLhs, _ConjRhs >
600 {
601 public:
602  typedef std::complex<RealScalar> Scalar;
603  typedef std::complex<RealScalar> LhsScalar;
604  typedef std::complex<RealScalar> RhsScalar;
605  typedef std::complex<RealScalar> ResScalar;
606 
607  enum {
608  ConjLhs = _ConjLhs,
609  ConjRhs = _ConjRhs,
610  Vectorizable = packet_traits<RealScalar>::Vectorizable
611  && packet_traits<Scalar>::Vectorizable,
612  RealPacketSize = Vectorizable ? packet_traits<RealScalar>::size : 1,
613  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
614  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
615  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
616 
617  // FIXME: should depend on NumberOfRegisters
618  nr = 4,
619  mr = ResPacketSize,
620 
621  LhsProgress = ResPacketSize,
622  RhsProgress = 1
623  };
624 
625  typedef typename packet_traits<RealScalar>::type RealPacket;
626  typedef typename packet_traits<Scalar>::type ScalarPacket;
627  typedef DoublePacket<RealPacket> DoublePacketType;
628 
629  typedef typename conditional<Vectorizable,RealPacket, Scalar>::type LhsPacket;
630  typedef typename conditional<Vectorizable,DoublePacketType,Scalar>::type RhsPacket;
631  typedef typename conditional<Vectorizable,ScalarPacket,Scalar>::type ResPacket;
632  typedef typename conditional<Vectorizable,DoublePacketType,Scalar>::type AccPacket;
633 
634  EIGEN_STRONG_INLINE void initAcc(Scalar& p) { p = Scalar(0); }
635 
636  EIGEN_STRONG_INLINE void initAcc(DoublePacketType& p)
637  {
638  p.first = pset1<RealPacket>(RealScalar(0));
639  p.second = pset1<RealPacket>(RealScalar(0));
640  }
641 
642  // Scalar path
643  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, ResPacket& dest) const
644  {
645  dest = pset1<ResPacket>(*b);
646  }
647 
648  // Vectorized path
649  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, DoublePacketType& dest) const
650  {
651  dest.first = pset1<RealPacket>(real(*b));
652  dest.second = pset1<RealPacket>(imag(*b));
653  }
654 
655  EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, ResPacket& dest) const
656  {
657  loadRhs(b,dest);
658  }
659  EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, DoublePacketType& dest) const
660  {
661  eigen_internal_assert(unpacket_traits<ScalarPacket>::size<=4);
662  loadRhs(b,dest);
663  }
664 
665  EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1, RhsPacket& b2, RhsPacket& b3)
666  {
667  // FIXME not sure that's the best way to implement it!
668  loadRhs(b+0, b0);
669  loadRhs(b+1, b1);
670  loadRhs(b+2, b2);
671  loadRhs(b+3, b3);
672  }
673 
674  // Vectorized path
675  EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, DoublePacketType& b0, DoublePacketType& b1)
676  {
677  // FIXME not sure that's the best way to implement it!
678  loadRhs(b+0, b0);
679  loadRhs(b+1, b1);
680  }
681 
682  // Scalar path
683  EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsScalar& b0, RhsScalar& b1)
684  {
685  // FIXME not sure that's the best way to implement it!
686  loadRhs(b+0, b0);
687  loadRhs(b+1, b1);
688  }
689 
690  // nothing special here
691  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
692  {
693  dest = pload<LhsPacket>((const typename unpacket_traits<LhsPacket>::type*)(a));
694  }
695 
696  EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacket& dest) const
697  {
698  dest = ploadu<LhsPacket>((const typename unpacket_traits<LhsPacket>::type*)(a));
699  }
700 
701  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, DoublePacketType& c, RhsPacket& /*tmp*/) const
702  {
703  c.first = padd(pmul(a,b.first), c.first);
704  c.second = padd(pmul(a,b.second),c.second);
705  }
706 
707  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, ResPacket& c, RhsPacket& /*tmp*/) const
708  {
709  c = cj.pmadd(a,b,c);
710  }
711 
712  EIGEN_STRONG_INLINE void acc(const Scalar& c, const Scalar& alpha, Scalar& r) const { r += alpha * c; }
713 
714  EIGEN_STRONG_INLINE void acc(const DoublePacketType& c, const ResPacket& alpha, ResPacket& r) const
715  {
716  // assemble c
717  ResPacket tmp;
718  if((!ConjLhs)&&(!ConjRhs))
719  {
720  tmp = pcplxflip(pconj(ResPacket(c.second)));
721  tmp = padd(ResPacket(c.first),tmp);
722  }
723  else if((!ConjLhs)&&(ConjRhs))
724  {
725  tmp = pconj(pcplxflip(ResPacket(c.second)));
726  tmp = padd(ResPacket(c.first),tmp);
727  }
728  else if((ConjLhs)&&(!ConjRhs))
729  {
730  tmp = pcplxflip(ResPacket(c.second));
731  tmp = padd(pconj(ResPacket(c.first)),tmp);
732  }
733  else if((ConjLhs)&&(ConjRhs))
734  {
735  tmp = pcplxflip(ResPacket(c.second));
736  tmp = psub(pconj(ResPacket(c.first)),tmp);
737  }
738 
739  r = pmadd(tmp,alpha,r);
740  }
741 
742 protected:
743  conj_helper<LhsScalar,RhsScalar,ConjLhs,ConjRhs> cj;
744 };
745 
746 template<typename RealScalar, bool _ConjRhs>
747 class gebp_traits<RealScalar, std::complex<RealScalar>, false, _ConjRhs >
748 {
749 public:
750  typedef std::complex<RealScalar> Scalar;
751  typedef RealScalar LhsScalar;
752  typedef Scalar RhsScalar;
753  typedef Scalar ResScalar;
754 
755  enum {
756  ConjLhs = false,
757  ConjRhs = _ConjRhs,
758  Vectorizable = packet_traits<RealScalar>::Vectorizable
759  && packet_traits<Scalar>::Vectorizable,
760  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
761  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
762  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1,
763 
764  NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
765  // FIXME: should depend on NumberOfRegisters
766  nr = 4,
767  mr = (EIGEN_PLAIN_ENUM_MIN(16,NumberOfRegisters)/2/nr)*ResPacketSize,
768 
769  LhsProgress = ResPacketSize,
770  RhsProgress = 1
771  };
772 
773  typedef typename packet_traits<LhsScalar>::type _LhsPacket;
774  typedef typename packet_traits<RhsScalar>::type _RhsPacket;
775  typedef typename packet_traits<ResScalar>::type _ResPacket;
776 
777  typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
778  typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
779  typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
780 
781  typedef ResPacket AccPacket;
782 
783  EIGEN_STRONG_INLINE void initAcc(AccPacket& p)
784  {
785  p = pset1<ResPacket>(ResScalar(0));
786  }
787 
788  EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
789  {
790  dest = pset1<RhsPacket>(*b);
791  }
792 
793  void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1, RhsPacket& b2, RhsPacket& b3)
794  {
795  pbroadcast4(b, b0, b1, b2, b3);
796  }
797 
798 // EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar* b, RhsPacket& b0, RhsPacket& b1)
799 // {
800 // // FIXME not sure that's the best way to implement it!
801 // b0 = pload1<RhsPacket>(b+0);
802 // b1 = pload1<RhsPacket>(b+1);
803 // }
804 
805  EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const
806  {
807  dest = ploaddup<LhsPacket>(a);
808  }
809 
810  EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const
811  {
812  eigen_internal_assert(unpacket_traits<RhsPacket>::size<=4);
813  loadRhs(b,dest);
814  }
815 
816  EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacket& dest) const
817  {
818  dest = ploaddup<LhsPacket>(a);
819  }
820 
821  EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp) const
822  {
823  madd_impl(a, b, c, tmp, typename conditional<Vectorizable,true_type,false_type>::type());
824  }
825 
826  EIGEN_STRONG_INLINE void madd_impl(const LhsPacket& a, const RhsPacket& b, AccPacket& c, RhsPacket& tmp, const true_type&) const
827  {
828 #ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
829  EIGEN_UNUSED_VARIABLE(tmp);
830  c.v = pmadd(a,b.v,c.v);
831 #else
832  tmp = b; tmp.v = pmul(a,tmp.v); c = padd(c,tmp);
833 #endif
834 
835  }
836 
837  EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/, const false_type&) const
838  {
839  c += a * b;
840  }
841 
842  EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const
843  {
844  r = cj.pmadd(alpha,c,r);
845  }
846 
847 protected:
848  conj_helper<ResPacket,ResPacket,false,ConjRhs> cj;
849 };
850 
851 /* optimized GEneral packed Block * packed Panel product kernel
852  *
853  * Mixing type logic: C += A * B
854  * | A | B | comments
855  * |real |cplx | no vectorization yet, would require to pack A with duplication
856  * |cplx |real | easy vectorization
857  */
858 template<typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
859 struct gebp_kernel
860 {
861  typedef gebp_traits<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> Traits;
862  typedef typename Traits::ResScalar ResScalar;
863  typedef typename Traits::LhsPacket LhsPacket;
864  typedef typename Traits::RhsPacket RhsPacket;
865  typedef typename Traits::ResPacket ResPacket;
866  typedef typename Traits::AccPacket AccPacket;
867 
868  typedef gebp_traits<RhsScalar,LhsScalar,ConjugateRhs,ConjugateLhs> SwappedTraits;
869  typedef typename SwappedTraits::ResScalar SResScalar;
870  typedef typename SwappedTraits::LhsPacket SLhsPacket;
871  typedef typename SwappedTraits::RhsPacket SRhsPacket;
872  typedef typename SwappedTraits::ResPacket SResPacket;
873  typedef typename SwappedTraits::AccPacket SAccPacket;
874 
875  typedef typename DataMapper::LinearMapper LinearMapper;
876 
877  enum {
878  Vectorizable = Traits::Vectorizable,
879  LhsProgress = Traits::LhsProgress,
880  RhsProgress = Traits::RhsProgress,
881  ResPacketSize = Traits::ResPacketSize
882  };
883 
884  EIGEN_DONT_INLINE
885  void operator()(const DataMapper& res, const LhsScalar* blockA, const RhsScalar* blockB,
886  Index rows, Index depth, Index cols, ResScalar alpha,
887  Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0);
888 };
889 
890 template<typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs>
891 EIGEN_DONT_INLINE
892 void gebp_kernel<LhsScalar,RhsScalar,Index,DataMapper,mr,nr,ConjugateLhs,ConjugateRhs>
893  ::operator()(const DataMapper& res, const LhsScalar* blockA, const RhsScalar* blockB,
894  Index rows, Index depth, Index cols, ResScalar alpha,
895  Index strideA, Index strideB, Index offsetA, Index offsetB)
896  {
897  Traits traits;
898  SwappedTraits straits;
899 
900  if(strideA==-1) strideA = depth;
901  if(strideB==-1) strideB = depth;
902  conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
903  Index packet_cols4 = nr>=4 ? (cols/4) * 4 : 0;
904  const Index peeled_mc3 = mr>=3*Traits::LhsProgress ? (rows/(3*LhsProgress))*(3*LhsProgress) : 0;
905  const Index peeled_mc2 = mr>=2*Traits::LhsProgress ? peeled_mc3+((rows-peeled_mc3)/(2*LhsProgress))*(2*LhsProgress) : 0;
906  const Index peeled_mc1 = mr>=1*Traits::LhsProgress ? (rows/(1*LhsProgress))*(1*LhsProgress) : 0;
907  enum { pk = 8 }; // NOTE Such a large peeling factor is important for large matrices (~ +5% when >1000 on Haswell)
908  const Index peeled_kc = depth & ~(pk-1);
909  const Index prefetch_res_offset = 32/sizeof(ResScalar);
910 // const Index depth2 = depth & ~1;
911 
912  //---------- Process 3 * LhsProgress rows at once ----------
913  // This corresponds to 3*LhsProgress x nr register blocks.
914  // Usually, make sense only with FMA
915  if(mr>=3*Traits::LhsProgress)
916  {
917  // Here, the general idea is to loop on each largest micro horizontal panel of the lhs (3*Traits::LhsProgress x depth)
918  // and on each largest micro vertical panel of the rhs (depth * nr).
919  // Blocking sizes, i.e., 'depth' has been computed so that the micro horizontal panel of the lhs fit in L1.
920  // However, if depth is too small, we can extend the number of rows of these horizontal panels.
921  // This actual number of rows is computed as follow:
922  const Index l1 = defaultL1CacheSize; // in Bytes, TODO, l1 should be passed to this function.
923  // The max(1, ...) here is needed because we may be using blocking params larger than what our known l1 cache size
924  // suggests we should be using: either because our known l1 cache size is inaccurate (e.g. on Android, we can only guess),
925  // or because we are testing specific blocking sizes.
926  const Index actual_panel_rows = (3*LhsProgress) * std::max<Index>(1,( (l1 - sizeof(ResScalar)*mr*nr - depth*nr*sizeof(RhsScalar)) / (depth * sizeof(LhsScalar) * 3*LhsProgress) ));
927  for(Index i1=0; i1<peeled_mc3; i1+=actual_panel_rows)
928  {
929  const Index actual_panel_end = (std::min)(i1+actual_panel_rows, peeled_mc3);
930  for(Index j2=0; j2<packet_cols4; j2+=nr)
931  {
932  for(Index i=i1; i<actual_panel_end; i+=3*LhsProgress)
933  {
934 
935  // We selected a 3*Traits::LhsProgress x nr micro block of res which is entirely
936  // stored into 3 x nr registers.
937 
938  const LhsScalar* blA = &blockA[i*strideA+offsetA*(3*LhsProgress)];
939  prefetch(&blA[0]);
940 
941  // gets res block as register
942  AccPacket C0, C1, C2, C3,
943  C4, C5, C6, C7,
944  C8, C9, C10, C11;
945  traits.initAcc(C0); traits.initAcc(C1); traits.initAcc(C2); traits.initAcc(C3);
946  traits.initAcc(C4); traits.initAcc(C5); traits.initAcc(C6); traits.initAcc(C7);
947  traits.initAcc(C8); traits.initAcc(C9); traits.initAcc(C10); traits.initAcc(C11);
948 
949  LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
950  LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
951  LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
952  LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
953 
954  r0.prefetch(0);
955  r1.prefetch(0);
956  r2.prefetch(0);
957  r3.prefetch(0);
958 
959  // performs "inner" products
960  const RhsScalar* blB = &blockB[j2*strideB+offsetB*nr];
961  prefetch(&blB[0]);
962  LhsPacket A0, A1;
963 
964  for(Index k=0; k<peeled_kc; k+=pk)
965  {
966  EIGEN_ASM_COMMENT("begin gebp micro kernel 3pX4");
967  RhsPacket B_0, T0;
968  LhsPacket A2;
969 
970 #define EIGEN_GEBP_ONESTEP(K) \
971  do { \
972  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 3pX4"); \
973  EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
974  internal::prefetch(blA+(3*K+16)*LhsProgress); \
975  if (EIGEN_ARCH_ARM) { internal::prefetch(blB+(4*K+16)*RhsProgress); } /* Bug 953 */ \
976  traits.loadLhs(&blA[(0+3*K)*LhsProgress], A0); \
977  traits.loadLhs(&blA[(1+3*K)*LhsProgress], A1); \
978  traits.loadLhs(&blA[(2+3*K)*LhsProgress], A2); \
979  traits.loadRhs(blB + (0+4*K)*Traits::RhsProgress, B_0); \
980  traits.madd(A0, B_0, C0, T0); \
981  traits.madd(A1, B_0, C4, T0); \
982  traits.madd(A2, B_0, C8, B_0); \
983  traits.loadRhs(blB + (1+4*K)*Traits::RhsProgress, B_0); \
984  traits.madd(A0, B_0, C1, T0); \
985  traits.madd(A1, B_0, C5, T0); \
986  traits.madd(A2, B_0, C9, B_0); \
987  traits.loadRhs(blB + (2+4*K)*Traits::RhsProgress, B_0); \
988  traits.madd(A0, B_0, C2, T0); \
989  traits.madd(A1, B_0, C6, T0); \
990  traits.madd(A2, B_0, C10, B_0); \
991  traits.loadRhs(blB + (3+4*K)*Traits::RhsProgress, B_0); \
992  traits.madd(A0, B_0, C3 , T0); \
993  traits.madd(A1, B_0, C7, T0); \
994  traits.madd(A2, B_0, C11, B_0); \
995  EIGEN_ASM_COMMENT("end step of gebp micro kernel 3pX4"); \
996  } while(false)
997 
998  internal::prefetch(blB);
999  EIGEN_GEBP_ONESTEP(0);
1000  EIGEN_GEBP_ONESTEP(1);
1001  EIGEN_GEBP_ONESTEP(2);
1002  EIGEN_GEBP_ONESTEP(3);
1003  EIGEN_GEBP_ONESTEP(4);
1004  EIGEN_GEBP_ONESTEP(5);
1005  EIGEN_GEBP_ONESTEP(6);
1006  EIGEN_GEBP_ONESTEP(7);
1007 
1008  blB += pk*4*RhsProgress;
1009  blA += pk*3*Traits::LhsProgress;
1010 
1011  EIGEN_ASM_COMMENT("end gebp micro kernel 3pX4");
1012  }
1013  // process remaining peeled loop
1014  for(Index k=peeled_kc; k<depth; k++)
1015  {
1016  RhsPacket B_0, T0;
1017  LhsPacket A2;
1018  EIGEN_GEBP_ONESTEP(0);
1019  blB += 4*RhsProgress;
1020  blA += 3*Traits::LhsProgress;
1021  }
1022 
1023 #undef EIGEN_GEBP_ONESTEP
1024 
1025  ResPacket R0, R1, R2;
1026  ResPacket alphav = pset1<ResPacket>(alpha);
1027 
1028  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1029  R1 = r0.loadPacket(1 * Traits::ResPacketSize);
1030  R2 = r0.loadPacket(2 * Traits::ResPacketSize);
1031  traits.acc(C0, alphav, R0);
1032  traits.acc(C4, alphav, R1);
1033  traits.acc(C8, alphav, R2);
1034  r0.storePacket(0 * Traits::ResPacketSize, R0);
1035  r0.storePacket(1 * Traits::ResPacketSize, R1);
1036  r0.storePacket(2 * Traits::ResPacketSize, R2);
1037 
1038  R0 = r1.loadPacket(0 * Traits::ResPacketSize);
1039  R1 = r1.loadPacket(1 * Traits::ResPacketSize);
1040  R2 = r1.loadPacket(2 * Traits::ResPacketSize);
1041  traits.acc(C1, alphav, R0);
1042  traits.acc(C5, alphav, R1);
1043  traits.acc(C9, alphav, R2);
1044  r1.storePacket(0 * Traits::ResPacketSize, R0);
1045  r1.storePacket(1 * Traits::ResPacketSize, R1);
1046  r1.storePacket(2 * Traits::ResPacketSize, R2);
1047 
1048  R0 = r2.loadPacket(0 * Traits::ResPacketSize);
1049  R1 = r2.loadPacket(1 * Traits::ResPacketSize);
1050  R2 = r2.loadPacket(2 * Traits::ResPacketSize);
1051  traits.acc(C2, alphav, R0);
1052  traits.acc(C6, alphav, R1);
1053  traits.acc(C10, alphav, R2);
1054  r2.storePacket(0 * Traits::ResPacketSize, R0);
1055  r2.storePacket(1 * Traits::ResPacketSize, R1);
1056  r2.storePacket(2 * Traits::ResPacketSize, R2);
1057 
1058  R0 = r3.loadPacket(0 * Traits::ResPacketSize);
1059  R1 = r3.loadPacket(1 * Traits::ResPacketSize);
1060  R2 = r3.loadPacket(2 * Traits::ResPacketSize);
1061  traits.acc(C3, alphav, R0);
1062  traits.acc(C7, alphav, R1);
1063  traits.acc(C11, alphav, R2);
1064  r3.storePacket(0 * Traits::ResPacketSize, R0);
1065  r3.storePacket(1 * Traits::ResPacketSize, R1);
1066  r3.storePacket(2 * Traits::ResPacketSize, R2);
1067  }
1068  }
1069 
1070  // Deal with remaining columns of the rhs
1071  for(Index j2=packet_cols4; j2<cols; j2++)
1072  {
1073  for(Index i=i1; i<actual_panel_end; i+=3*LhsProgress)
1074  {
1075  // One column at a time
1076  const LhsScalar* blA = &blockA[i*strideA+offsetA*(3*Traits::LhsProgress)];
1077  prefetch(&blA[0]);
1078 
1079  // gets res block as register
1080  AccPacket C0, C4, C8;
1081  traits.initAcc(C0);
1082  traits.initAcc(C4);
1083  traits.initAcc(C8);
1084 
1085  LinearMapper r0 = res.getLinearMapper(i, j2);
1086  r0.prefetch(0);
1087 
1088  // performs "inner" products
1089  const RhsScalar* blB = &blockB[j2*strideB+offsetB];
1090  LhsPacket A0, A1, A2;
1091 
1092  for(Index k=0; k<peeled_kc; k+=pk)
1093  {
1094  EIGEN_ASM_COMMENT("begin gebp micro kernel 3pX1");
1095  RhsPacket B_0;
1096 #define EIGEN_GEBGP_ONESTEP(K) \
1097  do { \
1098  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 3pX1"); \
1099  EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
1100  traits.loadLhs(&blA[(0+3*K)*LhsProgress], A0); \
1101  traits.loadLhs(&blA[(1+3*K)*LhsProgress], A1); \
1102  traits.loadLhs(&blA[(2+3*K)*LhsProgress], A2); \
1103  traits.loadRhs(&blB[(0+K)*RhsProgress], B_0); \
1104  traits.madd(A0, B_0, C0, B_0); \
1105  traits.madd(A1, B_0, C4, B_0); \
1106  traits.madd(A2, B_0, C8, B_0); \
1107  EIGEN_ASM_COMMENT("end step of gebp micro kernel 3pX1"); \
1108  } while(false)
1109 
1110  EIGEN_GEBGP_ONESTEP(0);
1111  EIGEN_GEBGP_ONESTEP(1);
1112  EIGEN_GEBGP_ONESTEP(2);
1113  EIGEN_GEBGP_ONESTEP(3);
1114  EIGEN_GEBGP_ONESTEP(4);
1115  EIGEN_GEBGP_ONESTEP(5);
1116  EIGEN_GEBGP_ONESTEP(6);
1117  EIGEN_GEBGP_ONESTEP(7);
1118 
1119  blB += pk*RhsProgress;
1120  blA += pk*3*Traits::LhsProgress;
1121 
1122  EIGEN_ASM_COMMENT("end gebp micro kernel 3pX1");
1123  }
1124 
1125  // process remaining peeled loop
1126  for(Index k=peeled_kc; k<depth; k++)
1127  {
1128  RhsPacket B_0;
1129  EIGEN_GEBGP_ONESTEP(0);
1130  blB += RhsProgress;
1131  blA += 3*Traits::LhsProgress;
1132  }
1133 #undef EIGEN_GEBGP_ONESTEP
1134  ResPacket R0, R1, R2;
1135  ResPacket alphav = pset1<ResPacket>(alpha);
1136 
1137  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1138  R1 = r0.loadPacket(1 * Traits::ResPacketSize);
1139  R2 = r0.loadPacket(2 * Traits::ResPacketSize);
1140  traits.acc(C0, alphav, R0);
1141  traits.acc(C4, alphav, R1);
1142  traits.acc(C8, alphav, R2);
1143  r0.storePacket(0 * Traits::ResPacketSize, R0);
1144  r0.storePacket(1 * Traits::ResPacketSize, R1);
1145  r0.storePacket(2 * Traits::ResPacketSize, R2);
1146  }
1147  }
1148  }
1149  }
1150 
1151  //---------- Process 2 * LhsProgress rows at once ----------
1152  if(mr>=2*Traits::LhsProgress)
1153  {
1154  const Index l1 = defaultL1CacheSize; // in Bytes, TODO, l1 should be passed to this function.
1155  // The max(1, ...) here is needed because we may be using blocking params larger than what our known l1 cache size
1156  // suggests we should be using: either because our known l1 cache size is inaccurate (e.g. on Android, we can only guess),
1157  // or because we are testing specific blocking sizes.
1158  Index actual_panel_rows = (2*LhsProgress) * std::max<Index>(1,( (l1 - sizeof(ResScalar)*mr*nr - depth*nr*sizeof(RhsScalar)) / (depth * sizeof(LhsScalar) * 2*LhsProgress) ));
1159 
1160  for(Index i1=peeled_mc3; i1<peeled_mc2; i1+=actual_panel_rows)
1161  {
1162  Index actual_panel_end = (std::min)(i1+actual_panel_rows, peeled_mc2);
1163  for(Index j2=0; j2<packet_cols4; j2+=nr)
1164  {
1165  for(Index i=i1; i<actual_panel_end; i+=2*LhsProgress)
1166  {
1167 
1168  // We selected a 2*Traits::LhsProgress x nr micro block of res which is entirely
1169  // stored into 2 x nr registers.
1170 
1171  const LhsScalar* blA = &blockA[i*strideA+offsetA*(2*Traits::LhsProgress)];
1172  prefetch(&blA[0]);
1173 
1174  // gets res block as register
1175  AccPacket C0, C1, C2, C3,
1176  C4, C5, C6, C7;
1177  traits.initAcc(C0); traits.initAcc(C1); traits.initAcc(C2); traits.initAcc(C3);
1178  traits.initAcc(C4); traits.initAcc(C5); traits.initAcc(C6); traits.initAcc(C7);
1179 
1180  LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
1181  LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
1182  LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
1183  LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
1184 
1185  r0.prefetch(prefetch_res_offset);
1186  r1.prefetch(prefetch_res_offset);
1187  r2.prefetch(prefetch_res_offset);
1188  r3.prefetch(prefetch_res_offset);
1189 
1190  // performs "inner" products
1191  const RhsScalar* blB = &blockB[j2*strideB+offsetB*nr];
1192  prefetch(&blB[0]);
1193  LhsPacket A0, A1;
1194 
1195  for(Index k=0; k<peeled_kc; k+=pk)
1196  {
1197  EIGEN_ASM_COMMENT("begin gebp micro kernel 2pX4");
1198  RhsPacket B_0, B1, B2, B3, T0;
1199 
1200  // NOTE: the begin/end asm comments below work around bug 935!
1201  // but they are not enough for gcc>=6 without FMA (bug 1637)
1202  #if EIGEN_GNUC_AT_LEAST(6,0) && defined(EIGEN_VECTORIZE_SSE)
1203  #define EIGEN_GEBP_2PX4_SPILLING_WORKAROUND __asm__ ("" : [a0] "+x,m" (A0),[a1] "+x,m" (A1));
1204  #else
1205  #define EIGEN_GEBP_2PX4_SPILLING_WORKAROUND
1206  #endif
1207  #define EIGEN_GEBGP_ONESTEP(K) \
1208  do { \
1209  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 2pX4"); \
1210  traits.loadLhs(&blA[(0+2*K)*LhsProgress], A0); \
1211  traits.loadLhs(&blA[(1+2*K)*LhsProgress], A1); \
1212  traits.broadcastRhs(&blB[(0+4*K)*RhsProgress], B_0, B1, B2, B3); \
1213  traits.madd(A0, B_0, C0, T0); \
1214  traits.madd(A1, B_0, C4, B_0); \
1215  traits.madd(A0, B1, C1, T0); \
1216  traits.madd(A1, B1, C5, B1); \
1217  traits.madd(A0, B2, C2, T0); \
1218  traits.madd(A1, B2, C6, B2); \
1219  traits.madd(A0, B3, C3, T0); \
1220  traits.madd(A1, B3, C7, B3); \
1221  EIGEN_GEBP_2PX4_SPILLING_WORKAROUND \
1222  EIGEN_ASM_COMMENT("end step of gebp micro kernel 2pX4"); \
1223  } while(false)
1224 
1225  internal::prefetch(blB+(48+0));
1226  EIGEN_GEBGP_ONESTEP(0);
1227  EIGEN_GEBGP_ONESTEP(1);
1228  EIGEN_GEBGP_ONESTEP(2);
1229  EIGEN_GEBGP_ONESTEP(3);
1230  internal::prefetch(blB+(48+16));
1231  EIGEN_GEBGP_ONESTEP(4);
1232  EIGEN_GEBGP_ONESTEP(5);
1233  EIGEN_GEBGP_ONESTEP(6);
1234  EIGEN_GEBGP_ONESTEP(7);
1235 
1236  blB += pk*4*RhsProgress;
1237  blA += pk*(2*Traits::LhsProgress);
1238 
1239  EIGEN_ASM_COMMENT("end gebp micro kernel 2pX4");
1240  }
1241  // process remaining peeled loop
1242  for(Index k=peeled_kc; k<depth; k++)
1243  {
1244  RhsPacket B_0, B1, B2, B3, T0;
1245  EIGEN_GEBGP_ONESTEP(0);
1246  blB += 4*RhsProgress;
1247  blA += 2*Traits::LhsProgress;
1248  }
1249 #undef EIGEN_GEBGP_ONESTEP
1250 
1251  ResPacket R0, R1, R2, R3;
1252  ResPacket alphav = pset1<ResPacket>(alpha);
1253 
1254  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1255  R1 = r0.loadPacket(1 * Traits::ResPacketSize);
1256  R2 = r1.loadPacket(0 * Traits::ResPacketSize);
1257  R3 = r1.loadPacket(1 * Traits::ResPacketSize);
1258  traits.acc(C0, alphav, R0);
1259  traits.acc(C4, alphav, R1);
1260  traits.acc(C1, alphav, R2);
1261  traits.acc(C5, alphav, R3);
1262  r0.storePacket(0 * Traits::ResPacketSize, R0);
1263  r0.storePacket(1 * Traits::ResPacketSize, R1);
1264  r1.storePacket(0 * Traits::ResPacketSize, R2);
1265  r1.storePacket(1 * Traits::ResPacketSize, R3);
1266 
1267  R0 = r2.loadPacket(0 * Traits::ResPacketSize);
1268  R1 = r2.loadPacket(1 * Traits::ResPacketSize);
1269  R2 = r3.loadPacket(0 * Traits::ResPacketSize);
1270  R3 = r3.loadPacket(1 * Traits::ResPacketSize);
1271  traits.acc(C2, alphav, R0);
1272  traits.acc(C6, alphav, R1);
1273  traits.acc(C3, alphav, R2);
1274  traits.acc(C7, alphav, R3);
1275  r2.storePacket(0 * Traits::ResPacketSize, R0);
1276  r2.storePacket(1 * Traits::ResPacketSize, R1);
1277  r3.storePacket(0 * Traits::ResPacketSize, R2);
1278  r3.storePacket(1 * Traits::ResPacketSize, R3);
1279  }
1280  }
1281 
1282  // Deal with remaining columns of the rhs
1283  for(Index j2=packet_cols4; j2<cols; j2++)
1284  {
1285  for(Index i=i1; i<actual_panel_end; i+=2*LhsProgress)
1286  {
1287  // One column at a time
1288  const LhsScalar* blA = &blockA[i*strideA+offsetA*(2*Traits::LhsProgress)];
1289  prefetch(&blA[0]);
1290 
1291  // gets res block as register
1292  AccPacket C0, C4;
1293  traits.initAcc(C0);
1294  traits.initAcc(C4);
1295 
1296  LinearMapper r0 = res.getLinearMapper(i, j2);
1297  r0.prefetch(prefetch_res_offset);
1298 
1299  // performs "inner" products
1300  const RhsScalar* blB = &blockB[j2*strideB+offsetB];
1301  LhsPacket A0, A1;
1302 
1303  for(Index k=0; k<peeled_kc; k+=pk)
1304  {
1305  EIGEN_ASM_COMMENT("begin gebp micro kernel 2pX1");
1306  RhsPacket B_0, B1;
1307 
1308 #define EIGEN_GEBGP_ONESTEP(K) \
1309  do { \
1310  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 2pX1"); \
1311  EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
1312  traits.loadLhs(&blA[(0+2*K)*LhsProgress], A0); \
1313  traits.loadLhs(&blA[(1+2*K)*LhsProgress], A1); \
1314  traits.loadRhs(&blB[(0+K)*RhsProgress], B_0); \
1315  traits.madd(A0, B_0, C0, B1); \
1316  traits.madd(A1, B_0, C4, B_0); \
1317  EIGEN_ASM_COMMENT("end step of gebp micro kernel 2pX1"); \
1318  } while(false)
1319 
1320  EIGEN_GEBGP_ONESTEP(0);
1321  EIGEN_GEBGP_ONESTEP(1);
1322  EIGEN_GEBGP_ONESTEP(2);
1323  EIGEN_GEBGP_ONESTEP(3);
1324  EIGEN_GEBGP_ONESTEP(4);
1325  EIGEN_GEBGP_ONESTEP(5);
1326  EIGEN_GEBGP_ONESTEP(6);
1327  EIGEN_GEBGP_ONESTEP(7);
1328 
1329  blB += pk*RhsProgress;
1330  blA += pk*2*Traits::LhsProgress;
1331 
1332  EIGEN_ASM_COMMENT("end gebp micro kernel 2pX1");
1333  }
1334 
1335  // process remaining peeled loop
1336  for(Index k=peeled_kc; k<depth; k++)
1337  {
1338  RhsPacket B_0, B1;
1339  EIGEN_GEBGP_ONESTEP(0);
1340  blB += RhsProgress;
1341  blA += 2*Traits::LhsProgress;
1342  }
1343 #undef EIGEN_GEBGP_ONESTEP
1344  ResPacket R0, R1;
1345  ResPacket alphav = pset1<ResPacket>(alpha);
1346 
1347  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1348  R1 = r0.loadPacket(1 * Traits::ResPacketSize);
1349  traits.acc(C0, alphav, R0);
1350  traits.acc(C4, alphav, R1);
1351  r0.storePacket(0 * Traits::ResPacketSize, R0);
1352  r0.storePacket(1 * Traits::ResPacketSize, R1);
1353  }
1354  }
1355  }
1356  }
1357  //---------- Process 1 * LhsProgress rows at once ----------
1358  if(mr>=1*Traits::LhsProgress)
1359  {
1360  // loops on each largest micro horizontal panel of lhs (1*LhsProgress x depth)
1361  for(Index i=peeled_mc2; i<peeled_mc1; i+=1*LhsProgress)
1362  {
1363  // loops on each largest micro vertical panel of rhs (depth * nr)
1364  for(Index j2=0; j2<packet_cols4; j2+=nr)
1365  {
1366  // We select a 1*Traits::LhsProgress x nr micro block of res which is entirely
1367  // stored into 1 x nr registers.
1368 
1369  const LhsScalar* blA = &blockA[i*strideA+offsetA*(1*Traits::LhsProgress)];
1370  prefetch(&blA[0]);
1371 
1372  // gets res block as register
1373  AccPacket C0, C1, C2, C3;
1374  traits.initAcc(C0);
1375  traits.initAcc(C1);
1376  traits.initAcc(C2);
1377  traits.initAcc(C3);
1378 
1379  LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
1380  LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
1381  LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
1382  LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
1383 
1384  r0.prefetch(prefetch_res_offset);
1385  r1.prefetch(prefetch_res_offset);
1386  r2.prefetch(prefetch_res_offset);
1387  r3.prefetch(prefetch_res_offset);
1388 
1389  // performs "inner" products
1390  const RhsScalar* blB = &blockB[j2*strideB+offsetB*nr];
1391  prefetch(&blB[0]);
1392  LhsPacket A0;
1393 
1394  for(Index k=0; k<peeled_kc; k+=pk)
1395  {
1396  EIGEN_ASM_COMMENT("begin gebp micro kernel 1pX4");
1397  RhsPacket B_0, B1, B2, B3;
1398 
1399 #define EIGEN_GEBGP_ONESTEP(K) \
1400  do { \
1401  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1pX4"); \
1402  EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
1403  traits.loadLhs(&blA[(0+1*K)*LhsProgress], A0); \
1404  traits.broadcastRhs(&blB[(0+4*K)*RhsProgress], B_0, B1, B2, B3); \
1405  traits.madd(A0, B_0, C0, B_0); \
1406  traits.madd(A0, B1, C1, B1); \
1407  traits.madd(A0, B2, C2, B2); \
1408  traits.madd(A0, B3, C3, B3); \
1409  EIGEN_ASM_COMMENT("end step of gebp micro kernel 1pX4"); \
1410  } while(false)
1411 
1412  internal::prefetch(blB+(48+0));
1413  EIGEN_GEBGP_ONESTEP(0);
1414  EIGEN_GEBGP_ONESTEP(1);
1415  EIGEN_GEBGP_ONESTEP(2);
1416  EIGEN_GEBGP_ONESTEP(3);
1417  internal::prefetch(blB+(48+16));
1418  EIGEN_GEBGP_ONESTEP(4);
1419  EIGEN_GEBGP_ONESTEP(5);
1420  EIGEN_GEBGP_ONESTEP(6);
1421  EIGEN_GEBGP_ONESTEP(7);
1422 
1423  blB += pk*4*RhsProgress;
1424  blA += pk*1*LhsProgress;
1425 
1426  EIGEN_ASM_COMMENT("end gebp micro kernel 1pX4");
1427  }
1428  // process remaining peeled loop
1429  for(Index k=peeled_kc; k<depth; k++)
1430  {
1431  RhsPacket B_0, B1, B2, B3;
1432  EIGEN_GEBGP_ONESTEP(0);
1433  blB += 4*RhsProgress;
1434  blA += 1*LhsProgress;
1435  }
1436 #undef EIGEN_GEBGP_ONESTEP
1437 
1438  ResPacket R0, R1;
1439  ResPacket alphav = pset1<ResPacket>(alpha);
1440 
1441  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1442  R1 = r1.loadPacket(0 * Traits::ResPacketSize);
1443  traits.acc(C0, alphav, R0);
1444  traits.acc(C1, alphav, R1);
1445  r0.storePacket(0 * Traits::ResPacketSize, R0);
1446  r1.storePacket(0 * Traits::ResPacketSize, R1);
1447 
1448  R0 = r2.loadPacket(0 * Traits::ResPacketSize);
1449  R1 = r3.loadPacket(0 * Traits::ResPacketSize);
1450  traits.acc(C2, alphav, R0);
1451  traits.acc(C3, alphav, R1);
1452  r2.storePacket(0 * Traits::ResPacketSize, R0);
1453  r3.storePacket(0 * Traits::ResPacketSize, R1);
1454  }
1455 
1456  // Deal with remaining columns of the rhs
1457  for(Index j2=packet_cols4; j2<cols; j2++)
1458  {
1459  // One column at a time
1460  const LhsScalar* blA = &blockA[i*strideA+offsetA*(1*Traits::LhsProgress)];
1461  prefetch(&blA[0]);
1462 
1463  // gets res block as register
1464  AccPacket C0;
1465  traits.initAcc(C0);
1466 
1467  LinearMapper r0 = res.getLinearMapper(i, j2);
1468 
1469  // performs "inner" products
1470  const RhsScalar* blB = &blockB[j2*strideB+offsetB];
1471  LhsPacket A0;
1472 
1473  for(Index k=0; k<peeled_kc; k+=pk)
1474  {
1475  EIGEN_ASM_COMMENT("begin gebp micro kernel 1pX1");
1476  RhsPacket B_0;
1477 
1478 #define EIGEN_GEBGP_ONESTEP(K) \
1479  do { \
1480  EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1pX1"); \
1481  EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
1482  traits.loadLhs(&blA[(0+1*K)*LhsProgress], A0); \
1483  traits.loadRhs(&blB[(0+K)*RhsProgress], B_0); \
1484  traits.madd(A0, B_0, C0, B_0); \
1485  EIGEN_ASM_COMMENT("end step of gebp micro kernel 1pX1"); \
1486  } while(false);
1487 
1488  EIGEN_GEBGP_ONESTEP(0);
1489  EIGEN_GEBGP_ONESTEP(1);
1490  EIGEN_GEBGP_ONESTEP(2);
1491  EIGEN_GEBGP_ONESTEP(3);
1492  EIGEN_GEBGP_ONESTEP(4);
1493  EIGEN_GEBGP_ONESTEP(5);
1494  EIGEN_GEBGP_ONESTEP(6);
1495  EIGEN_GEBGP_ONESTEP(7);
1496 
1497  blB += pk*RhsProgress;
1498  blA += pk*1*Traits::LhsProgress;
1499 
1500  EIGEN_ASM_COMMENT("end gebp micro kernel 1pX1");
1501  }
1502 
1503  // process remaining peeled loop
1504  for(Index k=peeled_kc; k<depth; k++)
1505  {
1506  RhsPacket B_0;
1507  EIGEN_GEBGP_ONESTEP(0);
1508  blB += RhsProgress;
1509  blA += 1*Traits::LhsProgress;
1510  }
1511 #undef EIGEN_GEBGP_ONESTEP
1512  ResPacket R0;
1513  ResPacket alphav = pset1<ResPacket>(alpha);
1514  R0 = r0.loadPacket(0 * Traits::ResPacketSize);
1515  traits.acc(C0, alphav, R0);
1516  r0.storePacket(0 * Traits::ResPacketSize, R0);
1517  }
1518  }
1519  }
1520  //---------- Process remaining rows, 1 at once ----------
1521  if(peeled_mc1<rows)
1522  {
1523  // loop on each panel of the rhs
1524  for(Index j2=0; j2<packet_cols4; j2+=nr)
1525  {
1526  // loop on each row of the lhs (1*LhsProgress x depth)
1527  for(Index i=peeled_mc1; i<rows; i+=1)
1528  {
1529  const LhsScalar* blA = &blockA[i*strideA+offsetA];
1530  prefetch(&blA[0]);
1531  const RhsScalar* blB = &blockB[j2*strideB+offsetB*nr];
1532 
1533  // The following piece of code wont work for 512 bit registers
1534  // Moreover, if LhsProgress==8 it assumes that there is a half packet of the same size
1535  // as nr (which is currently 4) for the return type.
1536  const int SResPacketHalfSize = unpacket_traits<typename unpacket_traits<SResPacket>::half>::size;
1537  if ((SwappedTraits::LhsProgress % 4) == 0 &&
1538  (SwappedTraits::LhsProgress <= 8) &&
1539  (SwappedTraits::LhsProgress!=8 || SResPacketHalfSize==nr))
1540  {
1541  SAccPacket C0, C1, C2, C3;
1542  straits.initAcc(C0);
1543  straits.initAcc(C1);
1544  straits.initAcc(C2);
1545  straits.initAcc(C3);
1546 
1547  const Index spk = (std::max)(1,SwappedTraits::LhsProgress/4);
1548  const Index endk = (depth/spk)*spk;
1549  const Index endk4 = (depth/(spk*4))*(spk*4);
1550 
1551  Index k=0;
1552  for(; k<endk4; k+=4*spk)
1553  {
1554  SLhsPacket A0,A1;
1555  SRhsPacket B_0,B_1;
1556 
1557  straits.loadLhsUnaligned(blB+0*SwappedTraits::LhsProgress, A0);
1558  straits.loadLhsUnaligned(blB+1*SwappedTraits::LhsProgress, A1);
1559 
1560  straits.loadRhsQuad(blA+0*spk, B_0);
1561  straits.loadRhsQuad(blA+1*spk, B_1);
1562  straits.madd(A0,B_0,C0,B_0);
1563  straits.madd(A1,B_1,C1,B_1);
1564 
1565  straits.loadLhsUnaligned(blB+2*SwappedTraits::LhsProgress, A0);
1566  straits.loadLhsUnaligned(blB+3*SwappedTraits::LhsProgress, A1);
1567  straits.loadRhsQuad(blA+2*spk, B_0);
1568  straits.loadRhsQuad(blA+3*spk, B_1);
1569  straits.madd(A0,B_0,C2,B_0);
1570  straits.madd(A1,B_1,C3,B_1);
1571 
1572  blB += 4*SwappedTraits::LhsProgress;
1573  blA += 4*spk;
1574  }
1575  C0 = padd(padd(C0,C1),padd(C2,C3));
1576  for(; k<endk; k+=spk)
1577  {
1578  SLhsPacket A0;
1579  SRhsPacket B_0;
1580 
1581  straits.loadLhsUnaligned(blB, A0);
1582  straits.loadRhsQuad(blA, B_0);
1583  straits.madd(A0,B_0,C0,B_0);
1584 
1585  blB += SwappedTraits::LhsProgress;
1586  blA += spk;
1587  }
1588  if(SwappedTraits::LhsProgress==8)
1589  {
1590  // Special case where we have to first reduce the accumulation register C0
1591  typedef typename conditional<SwappedTraits::LhsProgress>=8,typename unpacket_traits<SResPacket>::half,SResPacket>::type SResPacketHalf;
1592  typedef typename conditional<SwappedTraits::LhsProgress>=8,typename unpacket_traits<SLhsPacket>::half,SLhsPacket>::type SLhsPacketHalf;
1593  typedef typename conditional<SwappedTraits::LhsProgress>=8,typename unpacket_traits<SLhsPacket>::half,SRhsPacket>::type SRhsPacketHalf;
1594  typedef typename conditional<SwappedTraits::LhsProgress>=8,typename unpacket_traits<SAccPacket>::half,SAccPacket>::type SAccPacketHalf;
1595 
1596  SResPacketHalf R = res.template gatherPacket<SResPacketHalf>(i, j2);
1597  SResPacketHalf alphav = pset1<SResPacketHalf>(alpha);
1598 
1599  if(depth-endk>0)
1600  {
1601  // We have to handle the last row of the rhs which corresponds to a half-packet
1602  SLhsPacketHalf a0;
1603  SRhsPacketHalf b0;
1604  straits.loadLhsUnaligned(blB, a0);
1605  straits.loadRhs(blA, b0);
1606  SAccPacketHalf c0 = predux_downto4(C0);
1607  straits.madd(a0,b0,c0,b0);
1608  straits.acc(c0, alphav, R);
1609  }
1610  else
1611  {
1612  straits.acc(predux_downto4(C0), alphav, R);
1613  }
1614  res.scatterPacket(i, j2, R);
1615  }
1616  else
1617  {
1618  SResPacket R = res.template gatherPacket<SResPacket>(i, j2);
1619  SResPacket alphav = pset1<SResPacket>(alpha);
1620  straits.acc(C0, alphav, R);
1621  res.scatterPacket(i, j2, R);
1622  }
1623  }
1624  else // scalar path
1625  {
1626  // get a 1 x 4 res block as registers
1627  ResScalar C0(0), C1(0), C2(0), C3(0);
1628 
1629  for(Index k=0; k<depth; k++)
1630  {
1631  LhsScalar A0;
1632  RhsScalar B_0, B_1;
1633 
1634  A0 = blA[k];
1635 
1636  B_0 = blB[0];
1637  B_1 = blB[1];
1638  CJMADD(cj,A0,B_0,C0, B_0);
1639  CJMADD(cj,A0,B_1,C1, B_1);
1640 
1641  B_0 = blB[2];
1642  B_1 = blB[3];
1643  CJMADD(cj,A0,B_0,C2, B_0);
1644  CJMADD(cj,A0,B_1,C3, B_1);
1645 
1646  blB += 4;
1647  }
1648  res(i, j2 + 0) += alpha * C0;
1649  res(i, j2 + 1) += alpha * C1;
1650  res(i, j2 + 2) += alpha * C2;
1651  res(i, j2 + 3) += alpha * C3;
1652  }
1653  }
1654  }
1655  // remaining columns
1656  for(Index j2=packet_cols4; j2<cols; j2++)
1657  {
1658  // loop on each row of the lhs (1*LhsProgress x depth)
1659  for(Index i=peeled_mc1; i<rows; i+=1)
1660  {
1661  const LhsScalar* blA = &blockA[i*strideA+offsetA];
1662  prefetch(&blA[0]);
1663  // gets a 1 x 1 res block as registers
1664  ResScalar C0(0);
1665  const RhsScalar* blB = &blockB[j2*strideB+offsetB];
1666  for(Index k=0; k<depth; k++)
1667  {
1668  LhsScalar A0 = blA[k];
1669  RhsScalar B_0 = blB[k];
1670  CJMADD(cj, A0, B_0, C0, B_0);
1671  }
1672  res(i, j2) += alpha * C0;
1673  }
1674  }
1675  }
1676  }
1677 
1678 
1679 #undef CJMADD
1680 
1681 // pack a block of the lhs
1682 // The traversal is as follow (mr==4):
1683 // 0 4 8 12 ...
1684 // 1 5 9 13 ...
1685 // 2 6 10 14 ...
1686 // 3 7 11 15 ...
1687 //
1688 // 16 20 24 28 ...
1689 // 17 21 25 29 ...
1690 // 18 22 26 30 ...
1691 // 19 23 27 31 ...
1692 //
1693 // 32 33 34 35 ...
1694 // 36 36 38 39 ...
1695 template<typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, bool Conjugate, bool PanelMode>
1696 struct gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, ColMajor, Conjugate, PanelMode>
1697 {
1698  typedef typename DataMapper::LinearMapper LinearMapper;
1699  EIGEN_DONT_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
1700 };
1701 
1702 template<typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, bool Conjugate, bool PanelMode>
1703 EIGEN_DONT_INLINE void gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, ColMajor, Conjugate, PanelMode>
1704  ::operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
1705 {
1706  typedef typename packet_traits<Scalar>::type Packet;
1707  enum { PacketSize = packet_traits<Scalar>::size };
1708 
1709  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS");
1710  EIGEN_UNUSED_VARIABLE(stride);
1711  EIGEN_UNUSED_VARIABLE(offset);
1712  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1713  eigen_assert( ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) || (Pack1<=4) );
1714  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1715  Index count = 0;
1716 
1717  const Index peeled_mc3 = Pack1>=3*PacketSize ? (rows/(3*PacketSize))*(3*PacketSize) : 0;
1718  const Index peeled_mc2 = Pack1>=2*PacketSize ? peeled_mc3+((rows-peeled_mc3)/(2*PacketSize))*(2*PacketSize) : 0;
1719  const Index peeled_mc1 = Pack1>=1*PacketSize ? (rows/(1*PacketSize))*(1*PacketSize) : 0;
1720  const Index peeled_mc0 = Pack2>=1*PacketSize ? peeled_mc1
1721  : Pack2>1 ? (rows/Pack2)*Pack2 : 0;
1722 
1723  Index i=0;
1724 
1725  // Pack 3 packets
1726  if(Pack1>=3*PacketSize)
1727  {
1728  for(; i<peeled_mc3; i+=3*PacketSize)
1729  {
1730  if(PanelMode) count += (3*PacketSize) * offset;
1731 
1732  for(Index k=0; k<depth; k++)
1733  {
1734  Packet A, B, C;
1735  A = lhs.loadPacket(i+0*PacketSize, k);
1736  B = lhs.loadPacket(i+1*PacketSize, k);
1737  C = lhs.loadPacket(i+2*PacketSize, k);
1738  pstore(blockA+count, cj.pconj(A)); count+=PacketSize;
1739  pstore(blockA+count, cj.pconj(B)); count+=PacketSize;
1740  pstore(blockA+count, cj.pconj(C)); count+=PacketSize;
1741  }
1742  if(PanelMode) count += (3*PacketSize) * (stride-offset-depth);
1743  }
1744  }
1745  // Pack 2 packets
1746  if(Pack1>=2*PacketSize)
1747  {
1748  for(; i<peeled_mc2; i+=2*PacketSize)
1749  {
1750  if(PanelMode) count += (2*PacketSize) * offset;
1751 
1752  for(Index k=0; k<depth; k++)
1753  {
1754  Packet A, B;
1755  A = lhs.loadPacket(i+0*PacketSize, k);
1756  B = lhs.loadPacket(i+1*PacketSize, k);
1757  pstore(blockA+count, cj.pconj(A)); count+=PacketSize;
1758  pstore(blockA+count, cj.pconj(B)); count+=PacketSize;
1759  }
1760  if(PanelMode) count += (2*PacketSize) * (stride-offset-depth);
1761  }
1762  }
1763  // Pack 1 packets
1764  if(Pack1>=1*PacketSize)
1765  {
1766  for(; i<peeled_mc1; i+=1*PacketSize)
1767  {
1768  if(PanelMode) count += (1*PacketSize) * offset;
1769 
1770  for(Index k=0; k<depth; k++)
1771  {
1772  Packet A;
1773  A = lhs.loadPacket(i+0*PacketSize, k);
1774  pstore(blockA+count, cj.pconj(A));
1775  count+=PacketSize;
1776  }
1777  if(PanelMode) count += (1*PacketSize) * (stride-offset-depth);
1778  }
1779  }
1780  // Pack scalars
1781  if(Pack2<PacketSize && Pack2>1)
1782  {
1783  for(; i<peeled_mc0; i+=Pack2)
1784  {
1785  if(PanelMode) count += Pack2 * offset;
1786 
1787  for(Index k=0; k<depth; k++)
1788  for(Index w=0; w<Pack2; w++)
1789  blockA[count++] = cj(lhs(i+w, k));
1790 
1791  if(PanelMode) count += Pack2 * (stride-offset-depth);
1792  }
1793  }
1794  for(; i<rows; i++)
1795  {
1796  if(PanelMode) count += offset;
1797  for(Index k=0; k<depth; k++)
1798  blockA[count++] = cj(lhs(i, k));
1799  if(PanelMode) count += (stride-offset-depth);
1800  }
1801 }
1802 
1803 template<typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, bool Conjugate, bool PanelMode>
1804 struct gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, RowMajor, Conjugate, PanelMode>
1805 {
1806  typedef typename DataMapper::LinearMapper LinearMapper;
1807  EIGEN_DONT_INLINE void operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride=0, Index offset=0);
1808 };
1809 
1810 template<typename Scalar, typename Index, typename DataMapper, int Pack1, int Pack2, bool Conjugate, bool PanelMode>
1811 EIGEN_DONT_INLINE void gemm_pack_lhs<Scalar, Index, DataMapper, Pack1, Pack2, RowMajor, Conjugate, PanelMode>
1812  ::operator()(Scalar* blockA, const DataMapper& lhs, Index depth, Index rows, Index stride, Index offset)
1813 {
1814  typedef typename packet_traits<Scalar>::type Packet;
1815  enum { PacketSize = packet_traits<Scalar>::size };
1816 
1817  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS");
1818  EIGEN_UNUSED_VARIABLE(stride);
1819  EIGEN_UNUSED_VARIABLE(offset);
1820  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1821  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1822  Index count = 0;
1823 
1824 // const Index peeled_mc3 = Pack1>=3*PacketSize ? (rows/(3*PacketSize))*(3*PacketSize) : 0;
1825 // const Index peeled_mc2 = Pack1>=2*PacketSize ? peeled_mc3+((rows-peeled_mc3)/(2*PacketSize))*(2*PacketSize) : 0;
1826 // const Index peeled_mc1 = Pack1>=1*PacketSize ? (rows/(1*PacketSize))*(1*PacketSize) : 0;
1827 
1828  int pack = Pack1;
1829  Index i = 0;
1830  while(pack>0)
1831  {
1832  Index remaining_rows = rows-i;
1833  Index peeled_mc = i+(remaining_rows/pack)*pack;
1834  for(; i<peeled_mc; i+=pack)
1835  {
1836  if(PanelMode) count += pack * offset;
1837 
1838  const Index peeled_k = (depth/PacketSize)*PacketSize;
1839  Index k=0;
1840  if(pack>=PacketSize)
1841  {
1842  for(; k<peeled_k; k+=PacketSize)
1843  {
1844  for (Index m = 0; m < pack; m += PacketSize)
1845  {
1846  PacketBlock<Packet> kernel;
1847  for (int p = 0; p < PacketSize; ++p) kernel.packet[p] = lhs.loadPacket(i+p+m, k);
1848  ptranspose(kernel);
1849  for (int p = 0; p < PacketSize; ++p) pstore(blockA+count+m+(pack)*p, cj.pconj(kernel.packet[p]));
1850  }
1851  count += PacketSize*pack;
1852  }
1853  }
1854  for(; k<depth; k++)
1855  {
1856  Index w=0;
1857  for(; w<pack-3; w+=4)
1858  {
1859  Scalar a(cj(lhs(i+w+0, k))),
1860  b(cj(lhs(i+w+1, k))),
1861  c(cj(lhs(i+w+2, k))),
1862  d(cj(lhs(i+w+3, k)));
1863  blockA[count++] = a;
1864  blockA[count++] = b;
1865  blockA[count++] = c;
1866  blockA[count++] = d;
1867  }
1868  if(pack%4)
1869  for(;w<pack;++w)
1870  blockA[count++] = cj(lhs(i+w, k));
1871  }
1872 
1873  if(PanelMode) count += pack * (stride-offset-depth);
1874  }
1875 
1876  pack -= PacketSize;
1877  if(pack<Pack2 && (pack+PacketSize)!=Pack2)
1878  pack = Pack2;
1879  }
1880 
1881  for(; i<rows; i++)
1882  {
1883  if(PanelMode) count += offset;
1884  for(Index k=0; k<depth; k++)
1885  blockA[count++] = cj(lhs(i, k));
1886  if(PanelMode) count += (stride-offset-depth);
1887  }
1888 }
1889 
1890 // copy a complete panel of the rhs
1891 // this version is optimized for column major matrices
1892 // The traversal order is as follow: (nr==4):
1893 // 0 1 2 3 12 13 14 15 24 27
1894 // 4 5 6 7 16 17 18 19 25 28
1895 // 8 9 10 11 20 21 22 23 26 29
1896 // . . . . . . . . . .
1897 template<typename Scalar, typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
1898 struct gemm_pack_rhs<Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
1899 {
1900  typedef typename packet_traits<Scalar>::type Packet;
1901  typedef typename DataMapper::LinearMapper LinearMapper;
1902  enum { PacketSize = packet_traits<Scalar>::size };
1903  EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
1904 };
1905 
1906 template<typename Scalar, typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
1907 EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode>
1908  ::operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
1909 {
1910  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS COLMAJOR");
1911  EIGEN_UNUSED_VARIABLE(stride);
1912  EIGEN_UNUSED_VARIABLE(offset);
1913  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
1914  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
1915  Index packet_cols8 = nr>=8 ? (cols/8) * 8 : 0;
1916  Index packet_cols4 = nr>=4 ? (cols/4) * 4 : 0;
1917  Index count = 0;
1918  const Index peeled_k = (depth/PacketSize)*PacketSize;
1919 // if(nr>=8)
1920 // {
1921 // for(Index j2=0; j2<packet_cols8; j2+=8)
1922 // {
1923 // // skip what we have before
1924 // if(PanelMode) count += 8 * offset;
1925 // const Scalar* b0 = &rhs[(j2+0)*rhsStride];
1926 // const Scalar* b1 = &rhs[(j2+1)*rhsStride];
1927 // const Scalar* b2 = &rhs[(j2+2)*rhsStride];
1928 // const Scalar* b3 = &rhs[(j2+3)*rhsStride];
1929 // const Scalar* b4 = &rhs[(j2+4)*rhsStride];
1930 // const Scalar* b5 = &rhs[(j2+5)*rhsStride];
1931 // const Scalar* b6 = &rhs[(j2+6)*rhsStride];
1932 // const Scalar* b7 = &rhs[(j2+7)*rhsStride];
1933 // Index k=0;
1934 // if(PacketSize==8) // TODO enbale vectorized transposition for PacketSize==4
1935 // {
1936 // for(; k<peeled_k; k+=PacketSize) {
1937 // PacketBlock<Packet> kernel;
1938 // for (int p = 0; p < PacketSize; ++p) {
1939 // kernel.packet[p] = ploadu<Packet>(&rhs[(j2+p)*rhsStride+k]);
1940 // }
1941 // ptranspose(kernel);
1942 // for (int p = 0; p < PacketSize; ++p) {
1943 // pstoreu(blockB+count, cj.pconj(kernel.packet[p]));
1944 // count+=PacketSize;
1945 // }
1946 // }
1947 // }
1948 // for(; k<depth; k++)
1949 // {
1950 // blockB[count+0] = cj(b0[k]);
1951 // blockB[count+1] = cj(b1[k]);
1952 // blockB[count+2] = cj(b2[k]);
1953 // blockB[count+3] = cj(b3[k]);
1954 // blockB[count+4] = cj(b4[k]);
1955 // blockB[count+5] = cj(b5[k]);
1956 // blockB[count+6] = cj(b6[k]);
1957 // blockB[count+7] = cj(b7[k]);
1958 // count += 8;
1959 // }
1960 // // skip what we have after
1961 // if(PanelMode) count += 8 * (stride-offset-depth);
1962 // }
1963 // }
1964 
1965  if(nr>=4)
1966  {
1967  for(Index j2=packet_cols8; j2<packet_cols4; j2+=4)
1968  {
1969  // skip what we have before
1970  if(PanelMode) count += 4 * offset;
1971  const LinearMapper dm0 = rhs.getLinearMapper(0, j2 + 0);
1972  const LinearMapper dm1 = rhs.getLinearMapper(0, j2 + 1);
1973  const LinearMapper dm2 = rhs.getLinearMapper(0, j2 + 2);
1974  const LinearMapper dm3 = rhs.getLinearMapper(0, j2 + 3);
1975 
1976  Index k=0;
1977  if((PacketSize%4)==0) // TODO enable vectorized transposition for PacketSize==2 ??
1978  {
1979  for(; k<peeled_k; k+=PacketSize) {
1980  PacketBlock<Packet,(PacketSize%4)==0?4:PacketSize> kernel;
1981  kernel.packet[0] = dm0.loadPacket(k);
1982  kernel.packet[1%PacketSize] = dm1.loadPacket(k);
1983  kernel.packet[2%PacketSize] = dm2.loadPacket(k);
1984  kernel.packet[3%PacketSize] = dm3.loadPacket(k);
1985  ptranspose(kernel);
1986  pstoreu(blockB+count+0*PacketSize, cj.pconj(kernel.packet[0]));
1987  pstoreu(blockB+count+1*PacketSize, cj.pconj(kernel.packet[1%PacketSize]));
1988  pstoreu(blockB+count+2*PacketSize, cj.pconj(kernel.packet[2%PacketSize]));
1989  pstoreu(blockB+count+3*PacketSize, cj.pconj(kernel.packet[3%PacketSize]));
1990  count+=4*PacketSize;
1991  }
1992  }
1993  for(; k<depth; k++)
1994  {
1995  blockB[count+0] = cj(dm0(k));
1996  blockB[count+1] = cj(dm1(k));
1997  blockB[count+2] = cj(dm2(k));
1998  blockB[count+3] = cj(dm3(k));
1999  count += 4;
2000  }
2001  // skip what we have after
2002  if(PanelMode) count += 4 * (stride-offset-depth);
2003  }
2004  }
2005 
2006  // copy the remaining columns one at a time (nr==1)
2007  for(Index j2=packet_cols4; j2<cols; ++j2)
2008  {
2009  if(PanelMode) count += offset;
2010  const LinearMapper dm0 = rhs.getLinearMapper(0, j2);
2011  for(Index k=0; k<depth; k++)
2012  {
2013  blockB[count] = cj(dm0(k));
2014  count += 1;
2015  }
2016  if(PanelMode) count += (stride-offset-depth);
2017  }
2018 }
2019 
2020 // this version is optimized for row major matrices
2021 template<typename Scalar, typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
2022 struct gemm_pack_rhs<Scalar, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
2023 {
2024  typedef typename packet_traits<Scalar>::type Packet;
2025  typedef typename DataMapper::LinearMapper LinearMapper;
2026  enum { PacketSize = packet_traits<Scalar>::size };
2027  EIGEN_DONT_INLINE void operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride=0, Index offset=0);
2028 };
2029 
2030 template<typename Scalar, typename Index, typename DataMapper, int nr, bool Conjugate, bool PanelMode>
2031 EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode>
2032  ::operator()(Scalar* blockB, const DataMapper& rhs, Index depth, Index cols, Index stride, Index offset)
2033 {
2034  EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS ROWMAJOR");
2035  EIGEN_UNUSED_VARIABLE(stride);
2036  EIGEN_UNUSED_VARIABLE(offset);
2037  eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride));
2038  conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj;
2039  Index packet_cols8 = nr>=8 ? (cols/8) * 8 : 0;
2040  Index packet_cols4 = nr>=4 ? (cols/4) * 4 : 0;
2041  Index count = 0;
2042 
2043 // if(nr>=8)
2044 // {
2045 // for(Index j2=0; j2<packet_cols8; j2+=8)
2046 // {
2047 // // skip what we have before
2048 // if(PanelMode) count += 8 * offset;
2049 // for(Index k=0; k<depth; k++)
2050 // {
2051 // if (PacketSize==8) {
2052 // Packet A = ploadu<Packet>(&rhs[k*rhsStride + j2]);
2053 // pstoreu(blockB+count, cj.pconj(A));
2054 // } else if (PacketSize==4) {
2055 // Packet A = ploadu<Packet>(&rhs[k*rhsStride + j2]);
2056 // Packet B = ploadu<Packet>(&rhs[k*rhsStride + j2 + PacketSize]);
2057 // pstoreu(blockB+count, cj.pconj(A));
2058 // pstoreu(blockB+count+PacketSize, cj.pconj(B));
2059 // } else {
2060 // const Scalar* b0 = &rhs[k*rhsStride + j2];
2061 // blockB[count+0] = cj(b0[0]);
2062 // blockB[count+1] = cj(b0[1]);
2063 // blockB[count+2] = cj(b0[2]);
2064 // blockB[count+3] = cj(b0[3]);
2065 // blockB[count+4] = cj(b0[4]);
2066 // blockB[count+5] = cj(b0[5]);
2067 // blockB[count+6] = cj(b0[6]);
2068 // blockB[count+7] = cj(b0[7]);
2069 // }
2070 // count += 8;
2071 // }
2072 // // skip what we have after
2073 // if(PanelMode) count += 8 * (stride-offset-depth);
2074 // }
2075 // }
2076  if(nr>=4)
2077  {
2078  for(Index j2=packet_cols8; j2<packet_cols4; j2+=4)
2079  {
2080  // skip what we have before
2081  if(PanelMode) count += 4 * offset;
2082  for(Index k=0; k<depth; k++)
2083  {
2084  if (PacketSize==4) {
2085  Packet A = rhs.loadPacket(k, j2);
2086  pstoreu(blockB+count, cj.pconj(A));
2087  count += PacketSize;
2088  } else {
2089  const LinearMapper dm0 = rhs.getLinearMapper(k, j2);
2090  blockB[count+0] = cj(dm0(0));
2091  blockB[count+1] = cj(dm0(1));
2092  blockB[count+2] = cj(dm0(2));
2093  blockB[count+3] = cj(dm0(3));
2094  count += 4;
2095  }
2096  }
2097  // skip what we have after
2098  if(PanelMode) count += 4 * (stride-offset-depth);
2099  }
2100  }
2101  // copy the remaining columns one at a time (nr==1)
2102  for(Index j2=packet_cols4; j2<cols; ++j2)
2103  {
2104  if(PanelMode) count += offset;
2105  for(Index k=0; k<depth; k++)
2106  {
2107  blockB[count] = cj(rhs(k, j2));
2108  count += 1;
2109  }
2110  if(PanelMode) count += stride-offset-depth;
2111  }
2112 }
2113 
2114 } // end namespace internal
2115 
2118 inline std::ptrdiff_t l1CacheSize()
2119 {
2120  std::ptrdiff_t l1, l2, l3;
2121  internal::manage_caching_sizes(GetAction, &l1, &l2, &l3);
2122  return l1;
2123 }
2124 
2127 inline std::ptrdiff_t l2CacheSize()
2128 {
2129  std::ptrdiff_t l1, l2, l3;
2130  internal::manage_caching_sizes(GetAction, &l1, &l2, &l3);
2131  return l2;
2132 }
2133 
2137 inline std::ptrdiff_t l3CacheSize()
2138 {
2139  std::ptrdiff_t l1, l2, l3;
2140  internal::manage_caching_sizes(GetAction, &l1, &l2, &l3);
2141  return l3;
2142 }
2143 
2149 inline void setCpuCacheSizes(std::ptrdiff_t l1, std::ptrdiff_t l2, std::ptrdiff_t l3)
2150 {
2151  internal::manage_caching_sizes(SetAction, &l1, &l2, &l3);
2152 }
2153 
2154 } // end namespace Eigen
2155 
2156 #endif // EIGEN_GENERAL_BLOCK_PANEL_H
Eigen::internal::gebp_traits::LhsScalar
_LhsScalar LhsScalar
Definition: GeneralBlockPanelKernel.h:354
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::RhsPacket
conditional< Vectorizable, _RhsPacket, RhsScalar >::type RhsPacket
Definition: GeneralBlockPanelKernel.h:497
Eigen::internal::pbroadcast4
EIGEN_DEVICE_FUNC void pbroadcast4(const typename unpacket_traits< Packet >::type *a, Packet &a0, Packet &a1, Packet &a2, Packet &a3)
Definition: GenericPacketMath.h:256
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::madd_impl
EIGEN_STRONG_INLINE void madd_impl(const LhsPacket &a, const RhsPacket &b, AccPacket &c, RhsPacket &tmp, const true_type &) const
Definition: GeneralBlockPanelKernel.h:542
EIGEN_ALWAYS_INLINE
#define EIGEN_ALWAYS_INLINE
Definition: Macros.h:509
Eigen::internal::gebp_kernel::SResScalar
SwappedTraits::ResScalar SResScalar
Definition: GeneralBlockPanelKernel.h:869
Eigen::setCpuCacheSizes
void setCpuCacheSizes(std::ptrdiff_t l1, std::ptrdiff_t l2, std::ptrdiff_t l3)
Definition: GeneralBlockPanelKernel.h:2149
Eigen::internal::gemm_pack_rhs< Scalar, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode >::LinearMapper
DataMapper::LinearMapper LinearMapper
Definition: GeneralBlockPanelKernel.h:2025
Eigen::internal::gebp_traits::_RhsPacket
packet_traits< RhsScalar >::type _RhsPacket
Definition: GeneralBlockPanelKernel.h:387
EIGEN_GEBGP_ONESTEP
#define EIGEN_GEBGP_ONESTEP(K)
EIGEN_HAS_SINGLE_INSTRUCTION_MADD
#define EIGEN_HAS_SINGLE_INSTRUCTION_MADD
Definition: AltiVec/PacketMath.h:22
max
#define max(a, b)
Definition: datatypes.h:20
EIGEN_STRONG_INLINE
#define EIGEN_STRONG_INLINE
Definition: Macros.h:494
Eigen::internal::gebp_kernel::operator()
EIGEN_DONT_INLINE void operator()(const DataMapper &res, const LhsScalar *blockA, const RhsScalar *blockB, Index rows, Index depth, Index cols, ResScalar alpha, Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
Definition: GeneralBlockPanelKernel.h:893
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::AccPacket
conditional< Vectorizable, DoublePacketType, Scalar >::type AccPacket
Definition: GeneralBlockPanelKernel.h:632
Eigen::internal::gemm_pack_lhs< Scalar, Index, DataMapper, Pack1, Pack2, RowMajor, Conjugate, PanelMode >::LinearMapper
DataMapper::LinearMapper LinearMapper
Definition: GeneralBlockPanelKernel.h:1806
Eigen::internal::unpacket_traits::type
T type
Definition: XprHelper.h:160
Eigen::ColMajor
Definition: Constants.h:320
Eigen::internal::gemm_pack_lhs
Definition: BlasUtil.h:28
Eigen::internal::gebp_traits::ResScalar
ScalarBinaryOpTraits< LhsScalar, RhsScalar >::ReturnType ResScalar
Definition: GeneralBlockPanelKernel.h:356
Eigen::internal::gebp_madd_selector< CJ, T, T, T, T >::run
static EIGEN_ALWAYS_INLINE void run(const CJ &cj, T &a, T &b, T &c, T &t)
Definition: GeneralBlockPanelKernel.h:324
Eigen::internal::gebp_traits::LhsPacket
conditional< Vectorizable, _LhsPacket, LhsScalar >::type LhsPacket
Definition: GeneralBlockPanelKernel.h:390
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::loadRhs
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar *b, RhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:507
Eigen::internal::conj_helper
Definition: BlasUtil.h:61
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::Scalar
std::complex< RealScalar > Scalar
Definition: GeneralBlockPanelKernel.h:602
Eigen::internal::false_type
Definition: Meta.h:55
Eigen::internal::unpacket_traits
Definition: XprHelper.h:158
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadLhsUnaligned
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:696
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::loadRhsQuad
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar *b, RhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:512
Eigen::internal::gebp_madd
EIGEN_STRONG_INLINE void gebp_madd(const CJ &cj, A &a, B &b, C &c, T &t)
Definition: GeneralBlockPanelKernel.h:331
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::LhsPacket
conditional< Vectorizable, _LhsPacket, LhsScalar >::type LhsPacket
Definition: GeneralBlockPanelKernel.h:496
Eigen::internal::useSpecificBlockingSizes
bool useSpecificBlockingSizes(Index &k, Index &m, Index &n)
Definition: GeneralBlockPanelKernel.h:263
Eigen::internal::gebp_traits::_LhsPacket
packet_traits< LhsScalar >::type _LhsPacket
Definition: GeneralBlockPanelKernel.h:386
Eigen::internal::gebp_traits
Definition: GeneralBlockPanelKernel.h:19
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::internal::gebp_madd_selector::run
static EIGEN_ALWAYS_INLINE void run(const CJ &cj, A &a, B &b, C &c, T &)
Definition: GeneralBlockPanelKernel.h:317
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::loadRhsQuad
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar *b, RhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:810
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadLhs
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:691
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::initAcc
EIGEN_STRONG_INLINE void initAcc(Scalar &p)
Definition: GeneralBlockPanelKernel.h:634
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::RhsScalar
std::complex< RealScalar > RhsScalar
Definition: GeneralBlockPanelKernel.h:604
Eigen::internal::conj_helper::pmadd
EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar &x, const RhsScalar &y, const Scalar &c) const
Definition: BlasUtil.h:65
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadRhs
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar *b, DoublePacketType &dest) const
Definition: GeneralBlockPanelKernel.h:649
Eigen::internal::gebp_traits::AccPacket
ResPacket AccPacket
Definition: GeneralBlockPanelKernel.h:394
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadRhs
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar *b, ResPacket &dest) const
Definition: GeneralBlockPanelKernel.h:643
Eigen::internal::manage_caching_sizes_helper
std::ptrdiff_t manage_caching_sizes_helper(std::ptrdiff_t a, std::ptrdiff_t b)
Definition: GeneralBlockPanelKernel.h:23
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::ResPacket
conditional< Vectorizable, _ResPacket, ResScalar >::type ResPacket
Definition: GeneralBlockPanelKernel.h:498
Eigen::internal::gemm_pack_rhs< Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode >::LinearMapper
DataMapper::LinearMapper LinearMapper
Definition: GeneralBlockPanelKernel.h:1901
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::ResPacket
conditional< Vectorizable, _ResPacket, ResScalar >::type ResPacket
Definition: GeneralBlockPanelKernel.h:779
Eigen
Definition: common.h:73
Eigen::imag
DerType::Scalar imag(const AutoDiffScalar< DerType > &)
Definition: AutoDiffScalar.h:547
Eigen::internal::gebp_kernel::ResScalar
Traits::ResScalar ResScalar
Definition: GeneralBlockPanelKernel.h:862
utility::tuple::size
static constexpr size_t size(Tuple< Args... > &)
Provides access to the number of elements in a tuple as a compile-time constant expression.
Definition: TensorSyclTuple.h:131
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::broadcastRhs
EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar *b, RhsScalar &b0, RhsScalar &b1)
Definition: GeneralBlockPanelKernel.h:683
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::RhsPacket
conditional< Vectorizable, _RhsPacket, RhsScalar >::type RhsPacket
Definition: GeneralBlockPanelKernel.h:778
EIGEN_ASM_COMMENT
#define EIGEN_ASM_COMMENT(X)
Definition: Macros.h:624
c
RealScalar c
Definition: level1_cplx_impl.h:103
std
Definition: Half.h:150
Eigen::internal::gebp_traits::acc
EIGEN_STRONG_INLINE void acc(const AccPacket &c, const ResPacket &alpha, ResPacket &r) const
Definition: GeneralBlockPanelKernel.h:450
Eigen::internal::gebp_kernel::SAccPacket
SwappedTraits::AccPacket SAccPacket
Definition: GeneralBlockPanelKernel.h:873
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:823
Eigen::internal::manage_caching_sizes
void manage_caching_sizes(Action action, std::ptrdiff_t *l1, std::ptrdiff_t *l2, std::ptrdiff_t *l3)
Definition: GeneralBlockPanelKernel.h:55
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::madd
EIGEN_STRONG_INLINE void madd(const LhsPacket &a, const RhsPacket &b, ResPacket &c, RhsPacket &) const
Definition: GeneralBlockPanelKernel.h:707
Eigen::internal::packet_traits::type
T type
Definition: GenericPacketMath.h:98
Eigen::internal::gebp_kernel::RhsPacket
Traits::RhsPacket RhsPacket
Definition: GeneralBlockPanelKernel.h:864
Eigen::internal::defaultL3CacheSize
const std::ptrdiff_t defaultL3CacheSize
Definition: GeneralBlockPanelKernel.h:35
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::cj
conj_helper< LhsScalar, RhsScalar, ConjLhs, ConjRhs > cj
Definition: GeneralBlockPanelKernel.h:743
Eigen::internal::gemm_pack_rhs< Scalar, Index, DataMapper, nr, ColMajor, Conjugate, PanelMode >::Packet
packet_traits< Scalar >::type Packet
Definition: GeneralBlockPanelKernel.h:1900
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::loadRhs
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar *b, RhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:788
Eigen::internal::pcplxflip
EIGEN_STRONG_INLINE Packet2cf pcplxflip(const Packet2cf &x)
Definition: SSE/Complex.h:242
Eigen::internal::unpacket_traits< DoublePacket< Packet > >::half
DoublePacket< Packet > half
Definition: GeneralBlockPanelKernel.h:588
EIGEN_DONT_INLINE
#define EIGEN_DONT_INLINE
Definition: Macros.h:517
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::RealPacket
packet_traits< RealScalar >::type RealPacket
Definition: GeneralBlockPanelKernel.h:625
Eigen::l3CacheSize
std::ptrdiff_t l3CacheSize()
Definition: GeneralBlockPanelKernel.h:2137
Eigen::internal::gebp_traits::madd
EIGEN_STRONG_INLINE void madd(const LhsPacketType &a, const RhsPacketType &b, AccPacketType &c, AccPacketType &tmp) const
Definition: GeneralBlockPanelKernel.h:435
Eigen::internal::padd
EIGEN_DEVICE_FUNC Packet padd(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:151
CJMADD
#define CJMADD(CJ, A, B, C, T)
Definition: GeneralBlockPanelKernel.h:336
Eigen::internal::packet_traits
Definition: GenericPacketMath.h:96
Eigen::internal::gebp_kernel::SRhsPacket
SwappedTraits::RhsPacket SRhsPacket
Definition: GeneralBlockPanelKernel.h:871
Eigen::internal::gebp_kernel::SResPacket
SwappedTraits::ResPacket SResPacket
Definition: GeneralBlockPanelKernel.h:872
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::cj
conj_helper< ResPacket, ResPacket, false, ConjRhs > cj
Definition: GeneralBlockPanelKernel.h:848
A
MatrixType A(a, *n, *n, *lda)
Eigen::internal::gebp_kernel::Traits
gebp_traits< LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs > Traits
Definition: GeneralBlockPanelKernel.h:861
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: common.h:85
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::initAcc
EIGEN_STRONG_INLINE void initAcc(AccPacket &p)
Definition: GeneralBlockPanelKernel.h:783
Eigen::internal::gebp_traits::acc
EIGEN_STRONG_INLINE void acc(const ResPacketHalf &c, const ResPacketHalf &alpha, ResPacketHalf &r) const
Definition: GeneralBlockPanelKernel.h:456
Eigen::internal::gebp_traits::loadRhs
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar *b, RhsPacketType &dest) const
Definition: GeneralBlockPanelKernel.h:412
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::_RhsPacket
packet_traits< RhsScalar >::type _RhsPacket
Definition: GeneralBlockPanelKernel.h:493
Eigen::internal::pstoreu
EIGEN_DEVICE_FUNC void pstoreu(Scalar *to, const Packet &from)
Definition: GenericPacketMath.h:289
Eigen::internal::gebp_traits::broadcastRhs
EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar *b, RhsPacket &b0, RhsPacket &b1, RhsPacket &b2, RhsPacket &b3)
Definition: GeneralBlockPanelKernel.h:401
Eigen::l2CacheSize
std::ptrdiff_t l2CacheSize()
Definition: GeneralBlockPanelKernel.h:2127
Eigen::internal::conj_helper::pmul
EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar &x, const RhsScalar &y) const
Definition: BlasUtil.h:68
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::acc
EIGEN_STRONG_INLINE void acc(const DoublePacketType &c, const ResPacket &alpha, ResPacket &r) const
Definition: GeneralBlockPanelKernel.h:714
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::broadcastRhs
void broadcastRhs(const RhsScalar *b, RhsPacket &b0, RhsPacket &b1, RhsPacket &b2, RhsPacket &b3)
Definition: GeneralBlockPanelKernel.h:793
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::ResPacket
conditional< Vectorizable, ScalarPacket, Scalar >::type ResPacket
Definition: GeneralBlockPanelKernel.h:631
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::_RhsPacket
packet_traits< RhsScalar >::type _RhsPacket
Definition: GeneralBlockPanelKernel.h:774
Eigen::internal::ptranspose
EIGEN_STRONG_INLINE void ptranspose(PacketBlock< Packet2cf, 2 > &kernel)
Definition: AltiVec/Complex.h:242
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::_LhsPacket
packet_traits< LhsScalar >::type _LhsPacket
Definition: GeneralBlockPanelKernel.h:492
Eigen::internal::queryCacheSizes
void queryCacheSizes(int &l1, int &l2, int &l3)
Definition: Memory.h:936
Eigen::internal::conj_if
Definition: BlasUtil.h:43
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Eigen::internal::gebp_traits::ResPacket
conditional< Vectorizable, _ResPacket, ResScalar >::type ResPacket
Definition: GeneralBlockPanelKernel.h:392
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:579
Eigen::internal::evaluateProductBlockingSizesHeuristic
void evaluateProductBlockingSizesHeuristic(Index &k, Index &m, Index &n, Index num_threads=1)
Definition: GeneralBlockPanelKernel.h:93
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::broadcastRhs
EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar *b, DoublePacketType &b0, DoublePacketType &b1)
Definition: GeneralBlockPanelKernel.h:675
Eigen::numext::div_ceil
T div_ceil(const T &a, const T &b)
Definition: Meta.h:505
EIGEN_GEBP_ONESTEP
#define EIGEN_GEBP_ONESTEP(K)
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::madd_impl
EIGEN_STRONG_INLINE void madd_impl(const LhsScalar &a, const RhsScalar &b, ResScalar &c, RhsScalar &, const false_type &) const
Definition: GeneralBlockPanelKernel.h:837
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:125
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::acc
EIGEN_STRONG_INLINE void acc(const AccPacket &c, const ResPacket &alpha, ResPacket &r) const
Definition: GeneralBlockPanelKernel.h:842
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::LhsPacket
conditional< Vectorizable, RealPacket, Scalar >::type LhsPacket
Definition: GeneralBlockPanelKernel.h:629
Eigen::internal::gebp_traits::loadRhsQuad
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar *b, RhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:417
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::madd_impl
EIGEN_STRONG_INLINE void madd_impl(const LhsScalar &a, const RhsScalar &b, ResScalar &c, RhsScalar &, const false_type &) const
Definition: GeneralBlockPanelKernel.h:552
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::RhsScalar
RealScalar RhsScalar
Definition: GeneralBlockPanelKernel.h:468
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::loadLhs
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:805
Eigen::internal::gebp_traits::RhsScalar
_RhsScalar RhsScalar
Definition: GeneralBlockPanelKernel.h:355
EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS
#define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS
Definition: AltiVec/PacketMath.h:31
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::_ResPacket
packet_traits< ResScalar >::type _ResPacket
Definition: GeneralBlockPanelKernel.h:775
Eigen::internal::pconj
EIGEN_STRONG_INLINE Packet2cf pconj(const Packet2cf &a)
Definition: AltiVec/Complex.h:101
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::AccPacket
ResPacket AccPacket
Definition: GeneralBlockPanelKernel.h:781
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::AccPacket
ResPacket AccPacket
Definition: GeneralBlockPanelKernel.h:500
Scalar
SCALAR Scalar
Definition: common.h:84
Eigen::internal::gebp_madd_selector
Definition: GeneralBlockPanelKernel.h:316
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::loadLhsUnaligned
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:816
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::ScalarPacket
packet_traits< Scalar >::type ScalarPacket
Definition: GeneralBlockPanelKernel.h:626
b
Scalar * b
Definition: cholesky.cpp:56
Eigen::internal::pstore
EIGEN_DEVICE_FUNC void pstore(Scalar *to, const Packet &from)
Definition: GenericPacketMath.h:285
Eigen::internal::gebp_kernel
Definition: GeneralBlockPanelKernel.h:859
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadRhsQuad
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar *b, ResPacket &dest) const
Definition: GeneralBlockPanelKernel.h:655
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::madd
EIGEN_STRONG_INLINE void madd(const LhsPacket &a, const RhsPacket &b, AccPacket &c, RhsPacket &tmp) const
Definition: GeneralBlockPanelKernel.h:537
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::cj
conj_helper< ResPacket, ResPacket, ConjLhs, false > cj
Definition: GeneralBlockPanelKernel.h:563
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::broadcastRhs
EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar *b, RhsPacket &b0, RhsPacket &b1, RhsPacket &b2, RhsPacket &b3)
Definition: GeneralBlockPanelKernel.h:665
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::loadRhsQuad
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar *b, DoublePacketType &dest) const
Definition: GeneralBlockPanelKernel.h:659
Eigen::internal::gebp_kernel::AccPacket
Traits::AccPacket AccPacket
Definition: GeneralBlockPanelKernel.h:866
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::LhsScalar
RealScalar LhsScalar
Definition: GeneralBlockPanelKernel.h:751
Eigen::internal::defaultL2CacheSize
const std::ptrdiff_t defaultL2CacheSize
Definition: GeneralBlockPanelKernel.h:34
Eigen::internal::unpacket_traits::half
T half
Definition: XprHelper.h:161
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::ResScalar
Scalar ResScalar
Definition: GeneralBlockPanelKernel.h:753
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::LhsScalar
std::complex< RealScalar > LhsScalar
Definition: GeneralBlockPanelKernel.h:603
Eigen::internal::gebp_traits::loadLhsUnaligned
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar *a, LhsPacketType &dest) const
Definition: GeneralBlockPanelKernel.h:429
Eigen::internal::gebp_traits::initAcc
EIGEN_STRONG_INLINE void initAcc(AccPacket &p)
Definition: GeneralBlockPanelKernel.h:396
Eigen::internal::gebp_kernel::LhsPacket
Traits::LhsPacket LhsPacket
Definition: GeneralBlockPanelKernel.h:863
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::_ResPacket
packet_traits< ResScalar >::type _ResPacket
Definition: GeneralBlockPanelKernel.h:494
Eigen::internal::gebp_kernel::SLhsPacket
SwappedTraits::LhsPacket SLhsPacket
Definition: GeneralBlockPanelKernel.h:870
Eigen::internal::gebp_traits::RhsPacket
conditional< Vectorizable, _RhsPacket, RhsScalar >::type RhsPacket
Definition: GeneralBlockPanelKernel.h:391
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::acc
EIGEN_STRONG_INLINE void acc(const Scalar &c, const Scalar &alpha, Scalar &r) const
Definition: GeneralBlockPanelKernel.h:712
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::broadcastRhs
EIGEN_STRONG_INLINE void broadcastRhs(const RhsScalar *b, RhsPacket &b0, RhsPacket &b1, RhsPacket &b2, RhsPacket &b3)
Definition: GeneralBlockPanelKernel.h:527
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::LhsScalar
std::complex< RealScalar > LhsScalar
Definition: GeneralBlockPanelKernel.h:467
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::ResScalar
ScalarBinaryOpTraits< LhsScalar, RhsScalar >::ReturnType ResScalar
Definition: GeneralBlockPanelKernel.h:469
EIGEN_PLAIN_ENUM_MIN
#define EIGEN_PLAIN_ENUM_MIN(a, b)
Definition: Macros.h:875
Eigen::internal::defaultL1CacheSize
const std::ptrdiff_t defaultL1CacheSize
Definition: GeneralBlockPanelKernel.h:33
Eigen::internal::gebp_kernel::SwappedTraits
gebp_traits< RhsScalar, LhsScalar, ConjugateRhs, ConjugateLhs > SwappedTraits
Definition: GeneralBlockPanelKernel.h:868
Eigen::internal::computeProductBlockingSizes
void computeProductBlockingSizes(Index &k, Index &m, Index &n, Index num_threads=1)
Computes the blocking parameters for a m x k times k x n matrix product.
Definition: GeneralBlockPanelKernel.h:297
a
Scalar * a
Definition: cholesky.cpp:26
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::loadLhsUnaligned
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:522
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::acc
EIGEN_STRONG_INLINE void acc(const AccPacket &c, const ResPacket &alpha, ResPacket &r) const
Definition: GeneralBlockPanelKernel.h:557
Eigen::RowMajor
Definition: Constants.h:322
Eigen::internal::pmadd
EIGEN_STRONG_INLINE Packet4f pmadd(const Packet4f &a, const Packet4f &b, const Packet4f &c)
Definition: AltiVec/PacketMath.h:388
Eigen::internal::psub
EIGEN_DEVICE_FUNC Packet psub(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:156
Eigen::internal::gemm_pack_rhs< Scalar, Index, DataMapper, nr, RowMajor, Conjugate, PanelMode >::Packet
packet_traits< Scalar >::type Packet
Definition: GeneralBlockPanelKernel.h:2024
Eigen::internal::gebp_kernel::ResPacket
Traits::ResPacket ResPacket
Definition: GeneralBlockPanelKernel.h:865
Eigen::internal::gebp_kernel::LinearMapper
DataMapper::LinearMapper LinearMapper
Definition: GeneralBlockPanelKernel.h:875
Eigen::ScalarBinaryOpTraits
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition: XprHelper.h:766
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::ResScalar
std::complex< RealScalar > ResScalar
Definition: GeneralBlockPanelKernel.h:605
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::initAcc
EIGEN_STRONG_INLINE void initAcc(AccPacket &p)
Definition: GeneralBlockPanelKernel.h:502
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::DoublePacketType
DoublePacket< RealPacket > DoublePacketType
Definition: GeneralBlockPanelKernel.h:627
B
MatrixType B(b, *n, *nrhs, *ldb)
Eigen::internal::gebp_traits::_ResPacket
packet_traits< ResScalar >::type _ResPacket
Definition: GeneralBlockPanelKernel.h:388
eigen_internal_assert
#define eigen_internal_assert(x)
Definition: Macros.h:585
Eigen::internal::prefetch
EIGEN_DEVICE_FUNC void prefetch(const Scalar *addr)
Definition: GenericPacketMath.h:299
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::initAcc
EIGEN_STRONG_INLINE void initAcc(DoublePacketType &p)
Definition: GeneralBlockPanelKernel.h:636
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::madd_impl
EIGEN_STRONG_INLINE void madd_impl(const LhsPacket &a, const RhsPacket &b, AccPacket &c, RhsPacket &tmp, const true_type &) const
Definition: GeneralBlockPanelKernel.h:826
Eigen::internal::true_type
Definition: Meta.h:54
Eigen::real
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:545
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::RhsPacket
conditional< Vectorizable, DoublePacketType, Scalar >::type RhsPacket
Definition: GeneralBlockPanelKernel.h:630
n
PlainMatrixType mat * n
Definition: eigenvalues.cpp:41
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::madd
EIGEN_STRONG_INLINE void madd(const LhsPacket &a, const RhsPacket &b, AccPacket &c, RhsPacket &tmp) const
Definition: GeneralBlockPanelKernel.h:821
Eigen::internal::gebp_traits< std::complex< RealScalar >, std::complex< RealScalar >, _ConjLhs, _ConjRhs >::madd
EIGEN_STRONG_INLINE void madd(const LhsPacket &a, const RhsPacket &b, DoublePacketType &c, RhsPacket &) const
Definition: GeneralBlockPanelKernel.h:701
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::RhsScalar
Scalar RhsScalar
Definition: GeneralBlockPanelKernel.h:752
Eigen::internal::pmul
EIGEN_DEVICE_FUNC Packet pmul(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:170
Eigen::l1CacheSize
std::ptrdiff_t l1CacheSize()
Definition: GeneralBlockPanelKernel.h:2118
Eigen::internal::gemm_pack_lhs< Scalar, Index, DataMapper, Pack1, Pack2, ColMajor, Conjugate, PanelMode >::LinearMapper
DataMapper::LinearMapper LinearMapper
Definition: GeneralBlockPanelKernel.h:1698
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::LhsPacket
conditional< Vectorizable, _LhsPacket, LhsScalar >::type LhsPacket
Definition: GeneralBlockPanelKernel.h:777
Eigen::internal::predux_downto4
EIGEN_DEVICE_FUNC conditional<(unpacket_traits< Packet >::size%8)==0, typename unpacket_traits< Packet >::half, Packet >::type predux_downto4(const Packet &a)
Definition: GenericPacketMath.h:332
EIGEN_UNUSED_VARIABLE
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:618
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::Scalar
std::complex< RealScalar > Scalar
Definition: GeneralBlockPanelKernel.h:750
Eigen::internal::gebp_traits< RealScalar, std::complex< RealScalar >, false, _ConjRhs >::_LhsPacket
packet_traits< LhsScalar >::type _LhsPacket
Definition: GeneralBlockPanelKernel.h:773
Eigen::internal::gemm_pack_rhs
Definition: BlasUtil.h:25
Eigen::internal::gebp_traits< std::complex< RealScalar >, RealScalar, _ConjLhs, false >::loadLhs
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar *a, LhsPacket &dest) const
Definition: GeneralBlockPanelKernel.h:517
Eigen::internal::gebp_traits::loadLhs
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar *a, LhsPacketType &dest) const
Definition: GeneralBlockPanelKernel.h:423

control_box_rst
Author(s): Christoph Rösmann
autogenerated on Mon Feb 28 2022 22:06:52