Another set of UI messages fixes and tweaks! No functional changes.
[blender.git] / extern / Eigen3 / Eigen / src / Core / products / GeneralMatrixVector.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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_GENERAL_MATRIX_VECTOR_H
26 #define EIGEN_GENERAL_MATRIX_VECTOR_H
27
28 namespace internal {
29
30 /* Optimized col-major matrix * vector product:
31  * This algorithm processes 4 columns at onces that allows to both reduce
32  * the number of load/stores of the result by a factor 4 and to reduce
33  * the instruction dependency. Moreover, we know that all bands have the
34  * same alignment pattern.
35  *
36  * Mixing type logic: C += alpha * A * B
37  *  |  A  |  B  |alpha| comments
38  *  |real |cplx |cplx | no vectorization
39  *  |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization
40  *  |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp
41  *  |cplx |real |real | optimal case, vectorization possible via real-cplx mul
42  */
43 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs>
44 struct general_matrix_vector_product<Index,LhsScalar,ColMajor,ConjugateLhs,RhsScalar,ConjugateRhs>
45 {
46 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
47
48 enum {
49   Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable
50               && int(packet_traits<LhsScalar>::size)==int(packet_traits<RhsScalar>::size),
51   LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
52   RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
53   ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1
54 };
55
56 typedef typename packet_traits<LhsScalar>::type  _LhsPacket;
57 typedef typename packet_traits<RhsScalar>::type  _RhsPacket;
58 typedef typename packet_traits<ResScalar>::type  _ResPacket;
59
60 typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
61 typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
62 typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
63
64 EIGEN_DONT_INLINE static void run(
65   Index rows, Index cols,
66   const LhsScalar* lhs, Index lhsStride,
67   const RhsScalar* rhs, Index rhsIncr,
68   ResScalar* res, Index
69   #ifdef EIGEN_INTERNAL_DEBUGGING
70     resIncr
71   #endif
72   , RhsScalar alpha)
73 {
74   eigen_internal_assert(resIncr==1);
75   #ifdef _EIGEN_ACCUMULATE_PACKETS
76   #error _EIGEN_ACCUMULATE_PACKETS has already been defined
77   #endif
78   #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) \
79     pstore(&res[j], \
80       padd(pload<ResPacket>(&res[j]), \
81         padd( \
82           padd(pcj.pmul(EIGEN_CAT(ploa , A0)<LhsPacket>(&lhs0[j]),    ptmp0), \
83                   pcj.pmul(EIGEN_CAT(ploa , A13)<LhsPacket>(&lhs1[j]),   ptmp1)), \
84           padd(pcj.pmul(EIGEN_CAT(ploa , A2)<LhsPacket>(&lhs2[j]),    ptmp2), \
85                   pcj.pmul(EIGEN_CAT(ploa , A13)<LhsPacket>(&lhs3[j]),   ptmp3)) )))
86
87   conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
88   conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj;
89   if(ConjugateRhs)
90     alpha = conj(alpha);
91
92   enum { AllAligned = 0, EvenAligned, FirstAligned, NoneAligned };
93   const Index columnsAtOnce = 4;
94   const Index peels = 2;
95   const Index LhsPacketAlignedMask = LhsPacketSize-1;
96   const Index ResPacketAlignedMask = ResPacketSize-1;
97   const Index PeelAlignedMask = ResPacketSize*peels-1;
98   const Index size = rows;
99   
100   // How many coeffs of the result do we have to skip to be aligned.
101   // Here we assume data are at least aligned on the base scalar type.
102   Index alignedStart = first_aligned(res,size);
103   Index alignedSize = ResPacketSize>1 ? alignedStart + ((size-alignedStart) & ~ResPacketAlignedMask) : 0;
104   const Index peeledSize  = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
105
106   const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
107   Index alignmentPattern = alignmentStep==0 ? AllAligned
108                        : alignmentStep==(LhsPacketSize/2) ? EvenAligned
109                        : FirstAligned;
110
111   // we cannot assume the first element is aligned because of sub-matrices
112   const Index lhsAlignmentOffset = first_aligned(lhs,size);
113
114   // find how many columns do we have to skip to be aligned with the result (if possible)
115   Index skipColumns = 0;
116   // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
117   if( (size_t(lhs)%sizeof(LhsScalar)) || (size_t(res)%sizeof(ResScalar)) )
118   {
119     alignedSize = 0;
120     alignedStart = 0;
121   }
122   else if (LhsPacketSize>1)
123   {
124     eigen_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0 || size<LhsPacketSize);
125
126     while (skipColumns<LhsPacketSize &&
127           alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%LhsPacketSize))
128       ++skipColumns;
129     if (skipColumns==LhsPacketSize)
130     {
131       // nothing can be aligned, no need to skip any column
132       alignmentPattern = NoneAligned;
133       skipColumns = 0;
134     }
135     else
136     {
137       skipColumns = (std::min)(skipColumns,cols);
138       // note that the skiped columns are processed later.
139     }
140
141     eigen_internal_assert(  (alignmentPattern==NoneAligned)
142                       || (skipColumns + columnsAtOnce >= cols)
143                       || LhsPacketSize > size
144                       || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(LhsPacket))==0);
145   }
146   else if(Vectorizable)
147   {
148     alignedStart = 0;
149     alignedSize = size;
150     alignmentPattern = AllAligned;
151   }
152
153   Index offset1 = (FirstAligned && alignmentStep==1?3:1);
154   Index offset3 = (FirstAligned && alignmentStep==1?1:3);
155
156   Index columnBound = ((cols-skipColumns)/columnsAtOnce)*columnsAtOnce + skipColumns;
157   for (Index i=skipColumns; i<columnBound; i+=columnsAtOnce)
158   {
159     RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs[i*rhsIncr]),
160               ptmp1 = pset1<RhsPacket>(alpha*rhs[(i+offset1)*rhsIncr]),
161               ptmp2 = pset1<RhsPacket>(alpha*rhs[(i+2)*rhsIncr]),
162               ptmp3 = pset1<RhsPacket>(alpha*rhs[(i+offset3)*rhsIncr]);
163
164     // this helps a lot generating better binary code
165     const LhsScalar *lhs0 = lhs + i*lhsStride,     *lhs1 = lhs + (i+offset1)*lhsStride,
166                     *lhs2 = lhs + (i+2)*lhsStride, *lhs3 = lhs + (i+offset3)*lhsStride;
167
168     if (Vectorizable)
169     {
170       /* explicit vectorization */
171       // process initial unaligned coeffs
172       for (Index j=0; j<alignedStart; ++j)
173       {
174         res[j] = cj.pmadd(lhs0[j], pfirst(ptmp0), res[j]);
175         res[j] = cj.pmadd(lhs1[j], pfirst(ptmp1), res[j]);
176         res[j] = cj.pmadd(lhs2[j], pfirst(ptmp2), res[j]);
177         res[j] = cj.pmadd(lhs3[j], pfirst(ptmp3), res[j]);
178       }
179
180       if (alignedSize>alignedStart)
181       {
182         switch(alignmentPattern)
183         {
184           case AllAligned:
185             for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
186               _EIGEN_ACCUMULATE_PACKETS(d,d,d);
187             break;
188           case EvenAligned:
189             for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
190               _EIGEN_ACCUMULATE_PACKETS(d,du,d);
191             break;
192           case FirstAligned:
193             if(peels>1)
194             {
195               LhsPacket A00, A01, A02, A03, A10, A11, A12, A13;
196               ResPacket T0, T1;
197
198               A01 = pload<LhsPacket>(&lhs1[alignedStart-1]);
199               A02 = pload<LhsPacket>(&lhs2[alignedStart-2]);
200               A03 = pload<LhsPacket>(&lhs3[alignedStart-3]);
201
202               for (Index j = alignedStart; j<peeledSize; j+=peels*ResPacketSize)
203               {
204                 A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]);  palign<1>(A01,A11);
205                 A12 = pload<LhsPacket>(&lhs2[j-2+LhsPacketSize]);  palign<2>(A02,A12);
206                 A13 = pload<LhsPacket>(&lhs3[j-3+LhsPacketSize]);  palign<3>(A03,A13);
207
208                 A00 = pload<LhsPacket>(&lhs0[j]);
209                 A10 = pload<LhsPacket>(&lhs0[j+LhsPacketSize]);
210                 T0  = pcj.pmadd(A00, ptmp0, pload<ResPacket>(&res[j]));
211                 T1  = pcj.pmadd(A10, ptmp0, pload<ResPacket>(&res[j+ResPacketSize]));
212
213                 T0  = pcj.pmadd(A01, ptmp1, T0);
214                 A01 = pload<LhsPacket>(&lhs1[j-1+2*LhsPacketSize]);  palign<1>(A11,A01);
215                 T0  = pcj.pmadd(A02, ptmp2, T0);
216                 A02 = pload<LhsPacket>(&lhs2[j-2+2*LhsPacketSize]);  palign<2>(A12,A02);
217                 T0  = pcj.pmadd(A03, ptmp3, T0);
218                 pstore(&res[j],T0);
219                 A03 = pload<LhsPacket>(&lhs3[j-3+2*LhsPacketSize]);  palign<3>(A13,A03);
220                 T1  = pcj.pmadd(A11, ptmp1, T1);
221                 T1  = pcj.pmadd(A12, ptmp2, T1);
222                 T1  = pcj.pmadd(A13, ptmp3, T1);
223                 pstore(&res[j+ResPacketSize],T1);
224               }
225             }
226             for (Index j = peeledSize; j<alignedSize; j+=ResPacketSize)
227               _EIGEN_ACCUMULATE_PACKETS(d,du,du);
228             break;
229           default:
230             for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
231               _EIGEN_ACCUMULATE_PACKETS(du,du,du);
232             break;
233         }
234       }
235     } // end explicit vectorization
236
237     /* process remaining coeffs (or all if there is no explicit vectorization) */
238     for (Index j=alignedSize; j<size; ++j)
239     {
240       res[j] = cj.pmadd(lhs0[j], pfirst(ptmp0), res[j]);
241       res[j] = cj.pmadd(lhs1[j], pfirst(ptmp1), res[j]);
242       res[j] = cj.pmadd(lhs2[j], pfirst(ptmp2), res[j]);
243       res[j] = cj.pmadd(lhs3[j], pfirst(ptmp3), res[j]);
244     }
245   }
246
247   // process remaining first and last columns (at most columnsAtOnce-1)
248   Index end = cols;
249   Index start = columnBound;
250   do
251   {
252     for (Index k=start; k<end; ++k)
253     {
254       RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs[k*rhsIncr]);
255       const LhsScalar* lhs0 = lhs + k*lhsStride;
256
257       if (Vectorizable)
258       {
259         /* explicit vectorization */
260         // process first unaligned result's coeffs
261         for (Index j=0; j<alignedStart; ++j)
262           res[j] += cj.pmul(lhs0[j], pfirst(ptmp0));
263         // process aligned result's coeffs
264         if ((size_t(lhs0+alignedStart)%sizeof(LhsPacket))==0)
265           for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize)
266             pstore(&res[i], pcj.pmadd(ploadu<LhsPacket>(&lhs0[i]), ptmp0, pload<ResPacket>(&res[i])));
267         else
268           for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize)
269             pstore(&res[i], pcj.pmadd(ploadu<LhsPacket>(&lhs0[i]), ptmp0, pload<ResPacket>(&res[i])));
270       }
271
272       // process remaining scalars (or all if no explicit vectorization)
273       for (Index i=alignedSize; i<size; ++i)
274         res[i] += cj.pmul(lhs0[i], pfirst(ptmp0));
275     }
276     if (skipColumns)
277     {
278       start = 0;
279       end = skipColumns;
280       skipColumns = 0;
281     }
282     else
283       break;
284   } while(Vectorizable);
285   #undef _EIGEN_ACCUMULATE_PACKETS
286 }
287 };
288
289 /* Optimized row-major matrix * vector product:
290  * This algorithm processes 4 rows at onces that allows to both reduce
291  * the number of load/stores of the result by a factor 4 and to reduce
292  * the instruction dependency. Moreover, we know that all bands have the
293  * same alignment pattern.
294  *
295  * Mixing type logic:
296  *  - alpha is always a complex (or converted to a complex)
297  *  - no vectorization
298  */
299 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs>
300 struct general_matrix_vector_product<Index,LhsScalar,RowMajor,ConjugateLhs,RhsScalar,ConjugateRhs>
301 {
302 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
303
304 enum {
305   Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable
306               && int(packet_traits<LhsScalar>::size)==int(packet_traits<RhsScalar>::size),
307   LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
308   RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
309   ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1
310 };
311
312 typedef typename packet_traits<LhsScalar>::type  _LhsPacket;
313 typedef typename packet_traits<RhsScalar>::type  _RhsPacket;
314 typedef typename packet_traits<ResScalar>::type  _ResPacket;
315
316 typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
317 typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
318 typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
319   
320 EIGEN_DONT_INLINE static void run(
321   Index rows, Index cols,
322   const LhsScalar* lhs, Index lhsStride,
323   const RhsScalar* rhs, Index rhsIncr,
324   ResScalar* res, Index resIncr,
325   ResScalar alpha)
326 {
327   EIGEN_UNUSED_VARIABLE(rhsIncr);
328   eigen_internal_assert(rhsIncr==1);
329   #ifdef _EIGEN_ACCUMULATE_PACKETS
330   #error _EIGEN_ACCUMULATE_PACKETS has already been defined
331   #endif
332
333   #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) {\
334     RhsPacket b = pload<RhsPacket>(&rhs[j]); \
335     ptmp0 = pcj.pmadd(EIGEN_CAT(ploa,A0) <LhsPacket>(&lhs0[j]), b, ptmp0); \
336     ptmp1 = pcj.pmadd(EIGEN_CAT(ploa,A13)<LhsPacket>(&lhs1[j]), b, ptmp1); \
337     ptmp2 = pcj.pmadd(EIGEN_CAT(ploa,A2) <LhsPacket>(&lhs2[j]), b, ptmp2); \
338     ptmp3 = pcj.pmadd(EIGEN_CAT(ploa,A13)<LhsPacket>(&lhs3[j]), b, ptmp3); }
339
340   conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
341   conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj;
342
343   enum { AllAligned=0, EvenAligned=1, FirstAligned=2, NoneAligned=3 };
344   const Index rowsAtOnce = 4;
345   const Index peels = 2;
346   const Index RhsPacketAlignedMask = RhsPacketSize-1;
347   const Index LhsPacketAlignedMask = LhsPacketSize-1;
348   const Index PeelAlignedMask = RhsPacketSize*peels-1;
349   const Index depth = cols;
350
351   // How many coeffs of the result do we have to skip to be aligned.
352   // Here we assume data are at least aligned on the base scalar type
353   // if that's not the case then vectorization is discarded, see below.
354   Index alignedStart = first_aligned(rhs, depth);
355   Index alignedSize = RhsPacketSize>1 ? alignedStart + ((depth-alignedStart) & ~RhsPacketAlignedMask) : 0;
356   const Index peeledSize  = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart;
357
358   const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
359   Index alignmentPattern = alignmentStep==0 ? AllAligned
360                          : alignmentStep==(LhsPacketSize/2) ? EvenAligned
361                          : FirstAligned;
362
363   // we cannot assume the first element is aligned because of sub-matrices
364   const Index lhsAlignmentOffset = first_aligned(lhs,depth);
365
366   // find how many rows do we have to skip to be aligned with rhs (if possible)
367   Index skipRows = 0;
368   // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
369   if( (sizeof(LhsScalar)!=sizeof(RhsScalar)) || (size_t(lhs)%sizeof(LhsScalar)) || (size_t(rhs)%sizeof(RhsScalar)) )
370   {
371     alignedSize = 0;
372     alignedStart = 0;
373   }
374   else if (LhsPacketSize>1)
375   {
376     eigen_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0  || depth<LhsPacketSize);
377
378     while (skipRows<LhsPacketSize &&
379            alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%LhsPacketSize))
380       ++skipRows;
381     if (skipRows==LhsPacketSize)
382     {
383       // nothing can be aligned, no need to skip any column
384       alignmentPattern = NoneAligned;
385       skipRows = 0;
386     }
387     else
388     {
389       skipRows = (std::min)(skipRows,Index(rows));
390       // note that the skiped columns are processed later.
391     }
392     eigen_internal_assert(  alignmentPattern==NoneAligned
393                       || LhsPacketSize==1
394                       || (skipRows + rowsAtOnce >= rows)
395                       || LhsPacketSize > depth
396                       || (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(LhsPacket))==0);
397   }
398   else if(Vectorizable)
399   {
400     alignedStart = 0;
401     alignedSize = depth;
402     alignmentPattern = AllAligned;
403   }
404
405   Index offset1 = (FirstAligned && alignmentStep==1?3:1);
406   Index offset3 = (FirstAligned && alignmentStep==1?1:3);
407
408   Index rowBound = ((rows-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows;
409   for (Index i=skipRows; i<rowBound; i+=rowsAtOnce)
410   {
411     EIGEN_ALIGN16 ResScalar tmp0 = ResScalar(0);
412     ResScalar tmp1 = ResScalar(0), tmp2 = ResScalar(0), tmp3 = ResScalar(0);
413
414     // this helps the compiler generating good binary code
415     const LhsScalar *lhs0 = lhs + i*lhsStride,     *lhs1 = lhs + (i+offset1)*lhsStride,
416                     *lhs2 = lhs + (i+2)*lhsStride, *lhs3 = lhs + (i+offset3)*lhsStride;
417
418     if (Vectorizable)
419     {
420       /* explicit vectorization */
421       ResPacket ptmp0 = pset1<ResPacket>(ResScalar(0)), ptmp1 = pset1<ResPacket>(ResScalar(0)),
422                 ptmp2 = pset1<ResPacket>(ResScalar(0)), ptmp3 = pset1<ResPacket>(ResScalar(0));
423
424       // process initial unaligned coeffs
425       // FIXME this loop get vectorized by the compiler !
426       for (Index j=0; j<alignedStart; ++j)
427       {
428         RhsScalar b = rhs[j];
429         tmp0 += cj.pmul(lhs0[j],b); tmp1 += cj.pmul(lhs1[j],b);
430         tmp2 += cj.pmul(lhs2[j],b); tmp3 += cj.pmul(lhs3[j],b);
431       }
432
433       if (alignedSize>alignedStart)
434       {
435         switch(alignmentPattern)
436         {
437           case AllAligned:
438             for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
439               _EIGEN_ACCUMULATE_PACKETS(d,d,d);
440             break;
441           case EvenAligned:
442             for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
443               _EIGEN_ACCUMULATE_PACKETS(d,du,d);
444             break;
445           case FirstAligned:
446             if (peels>1)
447             {
448               /* Here we proccess 4 rows with with two peeled iterations to hide
449                * tghe overhead of unaligned loads. Moreover unaligned loads are handled
450                * using special shift/move operations between the two aligned packets
451                * overlaping the desired unaligned packet. This is *much* more efficient
452                * than basic unaligned loads.
453                */
454               LhsPacket A01, A02, A03, A11, A12, A13;
455               A01 = pload<LhsPacket>(&lhs1[alignedStart-1]);
456               A02 = pload<LhsPacket>(&lhs2[alignedStart-2]);
457               A03 = pload<LhsPacket>(&lhs3[alignedStart-3]);
458
459               for (Index j = alignedStart; j<peeledSize; j+=peels*RhsPacketSize)
460               {
461                 RhsPacket b = pload<RhsPacket>(&rhs[j]);
462                 A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]);  palign<1>(A01,A11);
463                 A12 = pload<LhsPacket>(&lhs2[j-2+LhsPacketSize]);  palign<2>(A02,A12);
464                 A13 = pload<LhsPacket>(&lhs3[j-3+LhsPacketSize]);  palign<3>(A03,A13);
465
466                 ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j]), b, ptmp0);
467                 ptmp1 = pcj.pmadd(A01, b, ptmp1);
468                 A01 = pload<LhsPacket>(&lhs1[j-1+2*LhsPacketSize]);  palign<1>(A11,A01);
469                 ptmp2 = pcj.pmadd(A02, b, ptmp2);
470                 A02 = pload<LhsPacket>(&lhs2[j-2+2*LhsPacketSize]);  palign<2>(A12,A02);
471                 ptmp3 = pcj.pmadd(A03, b, ptmp3);
472                 A03 = pload<LhsPacket>(&lhs3[j-3+2*LhsPacketSize]);  palign<3>(A13,A03);
473
474                 b = pload<RhsPacket>(&rhs[j+RhsPacketSize]);
475                 ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j+LhsPacketSize]), b, ptmp0);
476                 ptmp1 = pcj.pmadd(A11, b, ptmp1);
477                 ptmp2 = pcj.pmadd(A12, b, ptmp2);
478                 ptmp3 = pcj.pmadd(A13, b, ptmp3);
479               }
480             }
481             for (Index j = peeledSize; j<alignedSize; j+=RhsPacketSize)
482               _EIGEN_ACCUMULATE_PACKETS(d,du,du);
483             break;
484           default:
485             for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
486               _EIGEN_ACCUMULATE_PACKETS(du,du,du);
487             break;
488         }
489         tmp0 += predux(ptmp0);
490         tmp1 += predux(ptmp1);
491         tmp2 += predux(ptmp2);
492         tmp3 += predux(ptmp3);
493       }
494     } // end explicit vectorization
495
496     // process remaining coeffs (or all if no explicit vectorization)
497     // FIXME this loop get vectorized by the compiler !
498     for (Index j=alignedSize; j<depth; ++j)
499     {
500       RhsScalar b = rhs[j];
501       tmp0 += cj.pmul(lhs0[j],b); tmp1 += cj.pmul(lhs1[j],b);
502       tmp2 += cj.pmul(lhs2[j],b); tmp3 += cj.pmul(lhs3[j],b);
503     }
504     res[i*resIncr]            += alpha*tmp0;
505     res[(i+offset1)*resIncr]  += alpha*tmp1;
506     res[(i+2)*resIncr]        += alpha*tmp2;
507     res[(i+offset3)*resIncr]  += alpha*tmp3;
508   }
509
510   // process remaining first and last rows (at most columnsAtOnce-1)
511   Index end = rows;
512   Index start = rowBound;
513   do
514   {
515     for (Index i=start; i<end; ++i)
516     {
517       EIGEN_ALIGN16 ResScalar tmp0 = ResScalar(0);
518       ResPacket ptmp0 = pset1<ResPacket>(tmp0);
519       const LhsScalar* lhs0 = lhs + i*lhsStride;
520       // process first unaligned result's coeffs
521       // FIXME this loop get vectorized by the compiler !
522       for (Index j=0; j<alignedStart; ++j)
523         tmp0 += cj.pmul(lhs0[j], rhs[j]);
524
525       if (alignedSize>alignedStart)
526       {
527         // process aligned rhs coeffs
528         if ((size_t(lhs0+alignedStart)%sizeof(LhsPacket))==0)
529           for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize)
530             ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j]), pload<RhsPacket>(&rhs[j]), ptmp0);
531         else
532           for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize)
533             ptmp0 = pcj.pmadd(ploadu<LhsPacket>(&lhs0[j]), pload<RhsPacket>(&rhs[j]), ptmp0);
534         tmp0 += predux(ptmp0);
535       }
536
537       // process remaining scalars
538       // FIXME this loop get vectorized by the compiler !
539       for (Index j=alignedSize; j<depth; ++j)
540         tmp0 += cj.pmul(lhs0[j], rhs[j]);
541       res[i*resIncr] += alpha*tmp0;
542     }
543     if (skipRows)
544     {
545       start = 0;
546       end = skipRows;
547       skipRows = 0;
548     }
549     else
550       break;
551   } while(Vectorizable);
552
553   #undef _EIGEN_ACCUMULATE_PACKETS
554 }
555 };
556
557 } // end namespace internal
558
559 #endif // EIGEN_GENERAL_MATRIX_VECTOR_H