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