GeneralMatrixVector.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_MATRIX_VECTOR_H
11 #define EIGEN_GENERAL_MATRIX_VECTOR_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 /* Optimized col-major matrix * vector product:
18  * This algorithm processes 4 columns at onces that allows to both reduce
19  * the number of load/stores of the result by a factor 4 and to reduce
20  * the instruction dependency. Moreover, we know that all bands have the
21  * same alignment pattern.
22  *
23  * Mixing type logic: C += alpha * A * B
24  * | A | B |alpha| comments
25  * |real |cplx |cplx | no vectorization
26  * |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization
27  * |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp
28  * |cplx |real |real | optimal case, vectorization possible via real-cplx mul
29  *
30  * Accesses to the matrix coefficients follow the following logic:
31  *
32  * - if all columns have the same alignment then
33  * - if the columns have the same alignment as the result vector, then easy! (-> AllAligned case)
34  * - otherwise perform unaligned loads only (-> NoneAligned case)
35  * - otherwise
36  * - if even columns have the same alignment then
37  * // odd columns are guaranteed to have the same alignment too
38  * - if even or odd columns have the same alignment as the result, then
39  * // for a register size of 2 scalars, this is guarantee to be the case (e.g., SSE with double)
40  * - perform half aligned and half unaligned loads (-> EvenAligned case)
41  * - otherwise perform unaligned loads only (-> NoneAligned case)
42  * - otherwise, if the register size is 4 scalars (e.g., SSE with float) then
43  * - one over 4 consecutive columns is guaranteed to be aligned with the result vector,
44  * perform simple aligned loads for this column and aligned loads plus re-alignment for the other. (-> FirstAligned case)
45  * // this re-alignment is done by the palign function implemented for SSE in Eigen/src/Core/arch/SSE/PacketMath.h
46  * - otherwise,
47  * // if we get here, this means the register size is greater than 4 (e.g., AVX with floats),
48  * // we currently fall back to the NoneAligned case
49  *
50  * The same reasoning apply for the transposed case.
51  *
52  * The last case (PacketSize>4) could probably be improved by generalizing the FirstAligned case, but since we do not support AVX yet...
53  * One might also wonder why in the EvenAligned case we perform unaligned loads instead of using the aligned-loads plus re-alignment
54  * strategy as in the FirstAligned case. The reason is that we observed that unaligned loads on a 8 byte boundary are not too slow
55  * compared to unaligned loads on a 4 byte boundary.
56  *
57  */
58 template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
59 struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>
60 {
62 
63 enum {
66  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
67  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
68  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1
69 };
70 
74 
78 
79 EIGEN_DONT_INLINE static void run(
81  const LhsMapper& lhs,
82  const RhsMapper& rhs,
83  ResScalar* res, Index resIncr,
84  RhsScalar alpha);
85 };
86 
87 template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
90  const LhsMapper& lhs,
91  const RhsMapper& rhs,
92  ResScalar* res, Index resIncr,
93  RhsScalar alpha)
94 {
95  EIGEN_UNUSED_VARIABLE(resIncr);
96  eigen_internal_assert(resIncr==1);
97  #ifdef _EIGEN_ACCUMULATE_PACKETS
98  #error _EIGEN_ACCUMULATE_PACKETS has already been defined
99  #endif
100  #define _EIGEN_ACCUMULATE_PACKETS(Alignment0,Alignment13,Alignment2) \
101  pstore(&res[j], \
102  padd(pload<ResPacket>(&res[j]), \
103  padd( \
104  padd(pcj.pmul(lhs0.template load<LhsPacket, Alignment0>(j), ptmp0), \
105  pcj.pmul(lhs1.template load<LhsPacket, Alignment13>(j), ptmp1)), \
106  padd(pcj.pmul(lhs2.template load<LhsPacket, Alignment2>(j), ptmp2), \
107  pcj.pmul(lhs3.template load<LhsPacket, Alignment13>(j), ptmp3)) )))
108 
109  typedef typename LhsMapper::VectorMapper LhsScalars;
110 
113  if(ConjugateRhs)
114  alpha = numext::conj(alpha);
115 
116  enum { AllAligned = 0, EvenAligned, FirstAligned, NoneAligned };
117  const Index columnsAtOnce = 4;
118  const Index peels = 2;
119  const Index LhsPacketAlignedMask = LhsPacketSize-1;
120  const Index ResPacketAlignedMask = ResPacketSize-1;
121 // const Index PeelAlignedMask = ResPacketSize*peels-1;
122  const Index size = rows;
123 
124  const Index lhsStride = lhs.stride();
125 
126  // How many coeffs of the result do we have to skip to be aligned.
127  // Here we assume data are at least aligned on the base scalar type.
128  Index alignedStart = internal::first_default_aligned(res,size);
129  Index alignedSize = ResPacketSize>1 ? alignedStart + ((size-alignedStart) & ~ResPacketAlignedMask) : 0;
130  const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1;
131 
132  const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
133  Index alignmentPattern = alignmentStep==0 ? AllAligned
134  : alignmentStep==(LhsPacketSize/2) ? EvenAligned
135  : FirstAligned;
136 
137  // we cannot assume the first element is aligned because of sub-matrices
138  const Index lhsAlignmentOffset = lhs.firstAligned(size);
139 
140  // find how many columns do we have to skip to be aligned with the result (if possible)
141  Index skipColumns = 0;
142  // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
143  if( (lhsAlignmentOffset < 0) || (lhsAlignmentOffset == size) || (UIntPtr(res)%sizeof(ResScalar)) )
144  {
145  alignedSize = 0;
146  alignedStart = 0;
147  alignmentPattern = NoneAligned;
148  }
149  else if(LhsPacketSize > 4)
150  {
151  // TODO: extend the code to support aligned loads whenever possible when LhsPacketSize > 4.
152  // Currently, it seems to be better to perform unaligned loads anyway
153  alignmentPattern = NoneAligned;
154  }
155  else if (LhsPacketSize>1)
156  {
157  // eigen_internal_assert(size_t(firstLhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0 || size<LhsPacketSize);
158 
159  while (skipColumns<LhsPacketSize &&
160  alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%LhsPacketSize))
161  ++skipColumns;
162  if (skipColumns==LhsPacketSize)
163  {
164  // nothing can be aligned, no need to skip any column
165  alignmentPattern = NoneAligned;
166  skipColumns = 0;
167  }
168  else
169  {
170  skipColumns = (std::min)(skipColumns,cols);
171  // note that the skiped columns are processed later.
172  }
173 
174  /* eigen_internal_assert( (alignmentPattern==NoneAligned)
175  || (skipColumns + columnsAtOnce >= cols)
176  || LhsPacketSize > size
177  || (size_t(firstLhs+alignedStart+lhsStride*skipColumns)%sizeof(LhsPacket))==0);*/
178  }
179  else if(Vectorizable)
180  {
181  alignedStart = 0;
182  alignedSize = size;
183  alignmentPattern = AllAligned;
184  }
185 
186  const Index offset1 = (alignmentPattern==FirstAligned && alignmentStep==1)?3:1;
187  const Index offset3 = (alignmentPattern==FirstAligned && alignmentStep==1)?1:3;
188 
189  Index columnBound = ((cols-skipColumns)/columnsAtOnce)*columnsAtOnce + skipColumns;
190  for (Index i=skipColumns; i<columnBound; i+=columnsAtOnce)
191  {
192  RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs(i, 0)),
193  ptmp1 = pset1<RhsPacket>(alpha*rhs(i+offset1, 0)),
194  ptmp2 = pset1<RhsPacket>(alpha*rhs(i+2, 0)),
195  ptmp3 = pset1<RhsPacket>(alpha*rhs(i+offset3, 0));
196 
197  // this helps a lot generating better binary code
198  const LhsScalars lhs0 = lhs.getVectorMapper(0, i+0), lhs1 = lhs.getVectorMapper(0, i+offset1),
199  lhs2 = lhs.getVectorMapper(0, i+2), lhs3 = lhs.getVectorMapper(0, i+offset3);
200 
201  if (Vectorizable)
202  {
203  /* explicit vectorization */
204  // process initial unaligned coeffs
205  for (Index j=0; j<alignedStart; ++j)
206  {
207  res[j] = cj.pmadd(lhs0(j), pfirst(ptmp0), res[j]);
208  res[j] = cj.pmadd(lhs1(j), pfirst(ptmp1), res[j]);
209  res[j] = cj.pmadd(lhs2(j), pfirst(ptmp2), res[j]);
210  res[j] = cj.pmadd(lhs3(j), pfirst(ptmp3), res[j]);
211  }
212 
213  if (alignedSize>alignedStart)
214  {
215  switch(alignmentPattern)
216  {
217  case AllAligned:
218  for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
220  break;
221  case EvenAligned:
222  for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
224  break;
225  case FirstAligned:
226  {
227  Index j = alignedStart;
228  if(peels>1)
229  {
230  LhsPacket A00, A01, A02, A03, A10, A11, A12, A13;
231  ResPacket T0, T1;
232 
233  A01 = lhs1.template load<LhsPacket, Aligned>(alignedStart-1);
234  A02 = lhs2.template load<LhsPacket, Aligned>(alignedStart-2);
235  A03 = lhs3.template load<LhsPacket, Aligned>(alignedStart-3);
236 
237  for (; j<peeledSize; j+=peels*ResPacketSize)
238  {
239  A11 = lhs1.template load<LhsPacket, Aligned>(j-1+LhsPacketSize); palign<1>(A01,A11);
240  A12 = lhs2.template load<LhsPacket, Aligned>(j-2+LhsPacketSize); palign<2>(A02,A12);
241  A13 = lhs3.template load<LhsPacket, Aligned>(j-3+LhsPacketSize); palign<3>(A03,A13);
242 
243  A00 = lhs0.template load<LhsPacket, Aligned>(j);
244  A10 = lhs0.template load<LhsPacket, Aligned>(j+LhsPacketSize);
245  T0 = pcj.pmadd(A00, ptmp0, pload<ResPacket>(&res[j]));
246  T1 = pcj.pmadd(A10, ptmp0, pload<ResPacket>(&res[j+ResPacketSize]));
247 
248  T0 = pcj.pmadd(A01, ptmp1, T0);
249  A01 = lhs1.template load<LhsPacket, Aligned>(j-1+2*LhsPacketSize); palign<1>(A11,A01);
250  T0 = pcj.pmadd(A02, ptmp2, T0);
251  A02 = lhs2.template load<LhsPacket, Aligned>(j-2+2*LhsPacketSize); palign<2>(A12,A02);
252  T0 = pcj.pmadd(A03, ptmp3, T0);
253  pstore(&res[j],T0);
254  A03 = lhs3.template load<LhsPacket, Aligned>(j-3+2*LhsPacketSize); palign<3>(A13,A03);
255  T1 = pcj.pmadd(A11, ptmp1, T1);
256  T1 = pcj.pmadd(A12, ptmp2, T1);
257  T1 = pcj.pmadd(A13, ptmp3, T1);
258  pstore(&res[j+ResPacketSize],T1);
259  }
260  }
261  for (; j<alignedSize; j+=ResPacketSize)
263  break;
264  }
265  default:
266  for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize)
268  break;
269  }
270  }
271  } // end explicit vectorization
272 
273  /* process remaining coeffs (or all if there is no explicit vectorization) */
274  for (Index j=alignedSize; j<size; ++j)
275  {
276  res[j] = cj.pmadd(lhs0(j), pfirst(ptmp0), res[j]);
277  res[j] = cj.pmadd(lhs1(j), pfirst(ptmp1), res[j]);
278  res[j] = cj.pmadd(lhs2(j), pfirst(ptmp2), res[j]);
279  res[j] = cj.pmadd(lhs3(j), pfirst(ptmp3), res[j]);
280  }
281  }
282 
283  // process remaining first and last columns (at most columnsAtOnce-1)
284  Index end = cols;
285  Index start = columnBound;
286  do
287  {
288  for (Index k=start; k<end; ++k)
289  {
290  RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs(k, 0));
291  const LhsScalars lhs0 = lhs.getVectorMapper(0, k);
292 
293  if (Vectorizable)
294  {
295  /* explicit vectorization */
296  // process first unaligned result's coeffs
297  for (Index j=0; j<alignedStart; ++j)
298  res[j] += cj.pmul(lhs0(j), pfirst(ptmp0));
299  // process aligned result's coeffs
300  if (lhs0.template aligned<LhsPacket>(alignedStart))
301  for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize)
302  pstore(&res[i], pcj.pmadd(lhs0.template load<LhsPacket, Aligned>(i), ptmp0, pload<ResPacket>(&res[i])));
303  else
304  for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize)
305  pstore(&res[i], pcj.pmadd(lhs0.template load<LhsPacket, Unaligned>(i), ptmp0, pload<ResPacket>(&res[i])));
306  }
307 
308  // process remaining scalars (or all if no explicit vectorization)
309  for (Index i=alignedSize; i<size; ++i)
310  res[i] += cj.pmul(lhs0(i), pfirst(ptmp0));
311  }
312  if (skipColumns)
313  {
314  start = 0;
315  end = skipColumns;
316  skipColumns = 0;
317  }
318  else
319  break;
320  } while(Vectorizable);
321  #undef _EIGEN_ACCUMULATE_PACKETS
322 }
323 
324 /* Optimized row-major matrix * vector product:
325  * This algorithm processes 4 rows at onces that allows to both reduce
326  * the number of load/stores of the result by a factor 4 and to reduce
327  * the instruction dependency. Moreover, we know that all bands have the
328  * same alignment pattern.
329  *
330  * Mixing type logic:
331  * - alpha is always a complex (or converted to a complex)
332  * - no vectorization
333  */
334 template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
335 struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>
336 {
338 
339 enum {
342  LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1,
343  RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1,
344  ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1
345 };
346 
350 
354 
355 EIGEN_DONT_INLINE static void run(
356  Index rows, Index cols,
357  const LhsMapper& lhs,
358  const RhsMapper& rhs,
359  ResScalar* res, Index resIncr,
360  ResScalar alpha);
361 };
362 
363 template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
365  Index rows, Index cols,
366  const LhsMapper& lhs,
367  const RhsMapper& rhs,
368  ResScalar* res, Index resIncr,
370 {
371  eigen_internal_assert(rhs.stride()==1);
372 
373  #ifdef _EIGEN_ACCUMULATE_PACKETS
374  #error _EIGEN_ACCUMULATE_PACKETS has already been defined
375  #endif
376 
377  #define _EIGEN_ACCUMULATE_PACKETS(Alignment0,Alignment13,Alignment2) {\
378  RhsPacket b = rhs.getVectorMapper(j, 0).template load<RhsPacket, Aligned>(0); \
379  ptmp0 = pcj.pmadd(lhs0.template load<LhsPacket, Alignment0>(j), b, ptmp0); \
380  ptmp1 = pcj.pmadd(lhs1.template load<LhsPacket, Alignment13>(j), b, ptmp1); \
381  ptmp2 = pcj.pmadd(lhs2.template load<LhsPacket, Alignment2>(j), b, ptmp2); \
382  ptmp3 = pcj.pmadd(lhs3.template load<LhsPacket, Alignment13>(j), b, ptmp3); }
383 
386 
387  typedef typename LhsMapper::VectorMapper LhsScalars;
388 
389  enum { AllAligned=0, EvenAligned=1, FirstAligned=2, NoneAligned=3 };
390  const Index rowsAtOnce = 4;
391  const Index peels = 2;
392  const Index RhsPacketAlignedMask = RhsPacketSize-1;
393  const Index LhsPacketAlignedMask = LhsPacketSize-1;
394  const Index depth = cols;
395  const Index lhsStride = lhs.stride();
396 
397  // How many coeffs of the result do we have to skip to be aligned.
398  // Here we assume data are at least aligned on the base scalar type
399  // if that's not the case then vectorization is discarded, see below.
400  Index alignedStart = rhs.firstAligned(depth);
401  Index alignedSize = RhsPacketSize>1 ? alignedStart + ((depth-alignedStart) & ~RhsPacketAlignedMask) : 0;
402  const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1;
403 
404  const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0;
405  Index alignmentPattern = alignmentStep==0 ? AllAligned
406  : alignmentStep==(LhsPacketSize/2) ? EvenAligned
407  : FirstAligned;
408 
409  // we cannot assume the first element is aligned because of sub-matrices
410  const Index lhsAlignmentOffset = lhs.firstAligned(depth);
411  const Index rhsAlignmentOffset = rhs.firstAligned(rows);
412 
413  // find how many rows do we have to skip to be aligned with rhs (if possible)
414  Index skipRows = 0;
415  // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats)
416  if( (sizeof(LhsScalar)!=sizeof(RhsScalar)) ||
417  (lhsAlignmentOffset < 0) || (lhsAlignmentOffset == depth) ||
418  (rhsAlignmentOffset < 0) || (rhsAlignmentOffset == rows) )
419  {
420  alignedSize = 0;
421  alignedStart = 0;
422  alignmentPattern = NoneAligned;
423  }
424  else if(LhsPacketSize > 4)
425  {
426  // TODO: extend the code to support aligned loads whenever possible when LhsPacketSize > 4.
427  alignmentPattern = NoneAligned;
428  }
429  else if (LhsPacketSize>1)
430  {
431  // eigen_internal_assert(size_t(firstLhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0 || depth<LhsPacketSize);
432 
433  while (skipRows<LhsPacketSize &&
434  alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%LhsPacketSize))
435  ++skipRows;
436  if (skipRows==LhsPacketSize)
437  {
438  // nothing can be aligned, no need to skip any column
439  alignmentPattern = NoneAligned;
440  skipRows = 0;
441  }
442  else
443  {
444  skipRows = (std::min)(skipRows,Index(rows));
445  // note that the skiped columns are processed later.
446  }
447  /* eigen_internal_assert( alignmentPattern==NoneAligned
448  || LhsPacketSize==1
449  || (skipRows + rowsAtOnce >= rows)
450  || LhsPacketSize > depth
451  || (size_t(firstLhs+alignedStart+lhsStride*skipRows)%sizeof(LhsPacket))==0);*/
452  }
453  else if(Vectorizable)
454  {
455  alignedStart = 0;
456  alignedSize = depth;
457  alignmentPattern = AllAligned;
458  }
459 
460  const Index offset1 = (alignmentPattern==FirstAligned && alignmentStep==1)?3:1;
461  const Index offset3 = (alignmentPattern==FirstAligned && alignmentStep==1)?1:3;
462 
463  Index rowBound = ((rows-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows;
464  for (Index i=skipRows; i<rowBound; i+=rowsAtOnce)
465  {
466  // FIXME: what is the purpose of this EIGEN_ALIGN_DEFAULT ??
468  ResScalar tmp1 = ResScalar(0), tmp2 = ResScalar(0), tmp3 = ResScalar(0);
469 
470  // this helps the compiler generating good binary code
471  const LhsScalars lhs0 = lhs.getVectorMapper(i+0, 0), lhs1 = lhs.getVectorMapper(i+offset1, 0),
472  lhs2 = lhs.getVectorMapper(i+2, 0), lhs3 = lhs.getVectorMapper(i+offset3, 0);
473 
474  if (Vectorizable)
475  {
476  /* explicit vectorization */
477  ResPacket ptmp0 = pset1<ResPacket>(ResScalar(0)), ptmp1 = pset1<ResPacket>(ResScalar(0)),
478  ptmp2 = pset1<ResPacket>(ResScalar(0)), ptmp3 = pset1<ResPacket>(ResScalar(0));
479 
480  // process initial unaligned coeffs
481  // FIXME this loop get vectorized by the compiler !
482  for (Index j=0; j<alignedStart; ++j)
483  {
484  RhsScalar b = rhs(j, 0);
485  tmp0 += cj.pmul(lhs0(j),b); tmp1 += cj.pmul(lhs1(j),b);
486  tmp2 += cj.pmul(lhs2(j),b); tmp3 += cj.pmul(lhs3(j),b);
487  }
488 
489  if (alignedSize>alignedStart)
490  {
491  switch(alignmentPattern)
492  {
493  case AllAligned:
494  for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
496  break;
497  case EvenAligned:
498  for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
500  break;
501  case FirstAligned:
502  {
503  Index j = alignedStart;
504  if (peels>1)
505  {
506  /* Here we proccess 4 rows with with two peeled iterations to hide
507  * the overhead of unaligned loads. Moreover unaligned loads are handled
508  * using special shift/move operations between the two aligned packets
509  * overlaping the desired unaligned packet. This is *much* more efficient
510  * than basic unaligned loads.
511  */
512  LhsPacket A01, A02, A03, A11, A12, A13;
513  A01 = lhs1.template load<LhsPacket, Aligned>(alignedStart-1);
514  A02 = lhs2.template load<LhsPacket, Aligned>(alignedStart-2);
515  A03 = lhs3.template load<LhsPacket, Aligned>(alignedStart-3);
516 
517  for (; j<peeledSize; j+=peels*RhsPacketSize)
518  {
519  RhsPacket b = rhs.getVectorMapper(j, 0).template load<RhsPacket, Aligned>(0);
520  A11 = lhs1.template load<LhsPacket, Aligned>(j-1+LhsPacketSize); palign<1>(A01,A11);
521  A12 = lhs2.template load<LhsPacket, Aligned>(j-2+LhsPacketSize); palign<2>(A02,A12);
522  A13 = lhs3.template load<LhsPacket, Aligned>(j-3+LhsPacketSize); palign<3>(A03,A13);
523 
524  ptmp0 = pcj.pmadd(lhs0.template load<LhsPacket, Aligned>(j), b, ptmp0);
525  ptmp1 = pcj.pmadd(A01, b, ptmp1);
526  A01 = lhs1.template load<LhsPacket, Aligned>(j-1+2*LhsPacketSize); palign<1>(A11,A01);
527  ptmp2 = pcj.pmadd(A02, b, ptmp2);
528  A02 = lhs2.template load<LhsPacket, Aligned>(j-2+2*LhsPacketSize); palign<2>(A12,A02);
529  ptmp3 = pcj.pmadd(A03, b, ptmp3);
530  A03 = lhs3.template load<LhsPacket, Aligned>(j-3+2*LhsPacketSize); palign<3>(A13,A03);
531 
532  b = rhs.getVectorMapper(j+RhsPacketSize, 0).template load<RhsPacket, Aligned>(0);
533  ptmp0 = pcj.pmadd(lhs0.template load<LhsPacket, Aligned>(j+LhsPacketSize), b, ptmp0);
534  ptmp1 = pcj.pmadd(A11, b, ptmp1);
535  ptmp2 = pcj.pmadd(A12, b, ptmp2);
536  ptmp3 = pcj.pmadd(A13, b, ptmp3);
537  }
538  }
539  for (; j<alignedSize; j+=RhsPacketSize)
541  break;
542  }
543  default:
544  for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize)
546  break;
547  }
548  tmp0 += predux(ptmp0);
549  tmp1 += predux(ptmp1);
550  tmp2 += predux(ptmp2);
551  tmp3 += predux(ptmp3);
552  }
553  } // end explicit vectorization
554 
555  // process remaining coeffs (or all if no explicit vectorization)
556  // FIXME this loop get vectorized by the compiler !
557  for (Index j=alignedSize; j<depth; ++j)
558  {
559  RhsScalar b = rhs(j, 0);
560  tmp0 += cj.pmul(lhs0(j),b); tmp1 += cj.pmul(lhs1(j),b);
561  tmp2 += cj.pmul(lhs2(j),b); tmp3 += cj.pmul(lhs3(j),b);
562  }
563  res[i*resIncr] += alpha*tmp0;
564  res[(i+offset1)*resIncr] += alpha*tmp1;
565  res[(i+2)*resIncr] += alpha*tmp2;
566  res[(i+offset3)*resIncr] += alpha*tmp3;
567  }
568 
569  // process remaining first and last rows (at most columnsAtOnce-1)
570  Index end = rows;
571  Index start = rowBound;
572  do
573  {
574  for (Index i=start; i<end; ++i)
575  {
577  ResPacket ptmp0 = pset1<ResPacket>(tmp0);
578  const LhsScalars lhs0 = lhs.getVectorMapper(i, 0);
579  // process first unaligned result's coeffs
580  // FIXME this loop get vectorized by the compiler !
581  for (Index j=0; j<alignedStart; ++j)
582  tmp0 += cj.pmul(lhs0(j), rhs(j, 0));
583 
584  if (alignedSize>alignedStart)
585  {
586  // process aligned rhs coeffs
587  if (lhs0.template aligned<LhsPacket>(alignedStart))
588  for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize)
589  ptmp0 = pcj.pmadd(lhs0.template load<LhsPacket, Aligned>(j), rhs.getVectorMapper(j, 0).template load<RhsPacket, Aligned>(0), ptmp0);
590  else
591  for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize)
592  ptmp0 = pcj.pmadd(lhs0.template load<LhsPacket, Unaligned>(j), rhs.getVectorMapper(j, 0).template load<RhsPacket, Aligned>(0), ptmp0);
593  tmp0 += predux(ptmp0);
594  }
595 
596  // process remaining scalars
597  // FIXME this loop get vectorized by the compiler !
598  for (Index j=alignedSize; j<depth; ++j)
599  tmp0 += cj.pmul(lhs0(j), rhs(j, 0));
600  res[i*resIncr] += alpha*tmp0;
601  }
602  if (skipRows)
603  {
604  start = 0;
605  end = skipRows;
606  skipRows = 0;
607  }
608  else
609  break;
610  } while(Vectorizable);
611 
612  #undef _EIGEN_ACCUMULATE_PACKETS
613 }
614 
615 } // end namespace internal
616 
617 } // end namespace Eigen
618 
619 #endif // EIGEN_GENERAL_MATRIX_VECTOR_H
Scalar * b
Definition: benchVecAdd.cpp:17
return int(ret)+1
#define min(a, b)
Definition: datatypes.h:19
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
#define _EIGEN_ACCUMULATE_PACKETS(Alignment0, Alignment13, Alignment2)
EIGEN_DEVICE_FUNC unpacket_traits< Packet >::type predux(const Packet &a)
std::size_t UIntPtr
Definition: Meta.h:51
#define EIGEN_DONT_INLINE
Definition: Macros.h:517
static Index first_default_aligned(const DenseBase< Derived > &m)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
RealScalar alpha
EIGEN_DEVICE_FUNC unpacket_traits< Packet >::type pfirst(const Packet &a)
EIGEN_DEVICE_FUNC void pstore(Scalar *to, const Packet &from)
EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar &x, const RhsScalar &y, const Scalar &c) const
Definition: BlasUtil.h:65
#define EIGEN_ALIGN_MAX
Definition: Macros.h:757
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition: XprHelper.h:766
#define eigen_internal_assert(x)
Definition: Macros.h:585
Pose2 T1(M_PI/4.0, Point2(sqrt(0.5), sqrt(0.5)))
void run(Expr &expr, Dev &dev)
Definition: TensorSyclRun.h:33
std::ptrdiff_t j
EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar &x, const RhsScalar &y) const
Definition: BlasUtil.h:68
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:618
ScalarWithExceptions conj(const ScalarWithExceptions &x)
Definition: exceptions.cpp:74


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:06