Medial Code Documentation
Loading...
Searching...
No Matches
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// 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
13namespace 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 {
61 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
62
63 enum {
64 Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable
65 && int(packet_traits<LhsScalar>::size) == int(packet_traits<RhsScalar>::size),
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
71 typedef typename packet_traits<LhsScalar>::type _LhsPacket;
72 typedef typename packet_traits<RhsScalar>::type _RhsPacket;
73 typedef typename packet_traits<ResScalar>::type _ResPacket;
74
75 typedef typename conditional<Vectorizable, _LhsPacket, LhsScalar>::type LhsPacket;
76 typedef typename conditional<Vectorizable, _RhsPacket, RhsScalar>::type RhsPacket;
77 typedef typename conditional<Vectorizable, _ResPacket, ResScalar>::type ResPacket;
78
79 EIGEN_DONT_INLINE static void run(
80 Index rows, Index cols,
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>
88 EIGEN_DONT_INLINE void general_matrix_vector_product<Index, LhsScalar, LhsMapper, ColMajor, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::run(
89 Index rows, Index cols,
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
111 conj_helper<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs> cj;
112 conj_helper<LhsPacket, RhsPacket, ConjugateLhs, ConjugateRhs> pcj;
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) || (size_t(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 = (FirstAligned != 0 && alignmentStep == 1 ? 3 : 1);
187 const Index offset3 = (FirstAligned != 0 && 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)
219 _EIGEN_ACCUMULATE_PACKETS(Aligned, Aligned, Aligned);
220 break;
221 case EvenAligned:
222 for (Index j = alignedStart; j < alignedSize; j += ResPacketSize)
223 _EIGEN_ACCUMULATE_PACKETS(Aligned, Unaligned, Aligned);
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)
262 _EIGEN_ACCUMULATE_PACKETS(Aligned, Unaligned, Unaligned);
263 break;
264 }
265 default:
266 for (Index j = alignedStart; j < alignedSize; j += ResPacketSize)
267 _EIGEN_ACCUMULATE_PACKETS(Unaligned, Unaligned, Unaligned);
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 {
337 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
338
339 enum {
340 Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable
341 && int(packet_traits<LhsScalar>::size) == int(packet_traits<RhsScalar>::size),
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
347 typedef typename packet_traits<LhsScalar>::type _LhsPacket;
348 typedef typename packet_traits<RhsScalar>::type _RhsPacket;
349 typedef typename packet_traits<ResScalar>::type _ResPacket;
350
351 typedef typename conditional<Vectorizable, _LhsPacket, LhsScalar>::type LhsPacket;
352 typedef typename conditional<Vectorizable, _RhsPacket, RhsScalar>::type RhsPacket;
353 typedef typename conditional<Vectorizable, _ResPacket, ResScalar>::type ResPacket;
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>
364 EIGEN_DONT_INLINE void general_matrix_vector_product<Index, LhsScalar, LhsMapper, RowMajor, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::run(
365 Index rows, Index cols,
366 const LhsMapper& lhs,
367 const RhsMapper& rhs,
368 ResScalar* res, Index resIncr,
369 ResScalar alpha)
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
384 conj_helper<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs> cj;
385 conj_helper<LhsPacket, RhsPacket, ConjugateLhs, ConjugateRhs> pcj;
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 = (FirstAligned != 0 && alignmentStep == 1 ? 3 : 1);
461 const Index offset3 = (FirstAligned != 0 && 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 ??
467 EIGEN_ALIGN_MAX ResScalar tmp0 = ResScalar(0);
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)
495 _EIGEN_ACCUMULATE_PACKETS(Aligned, Aligned, Aligned);
496 break;
497 case EvenAligned:
498 for (Index j = alignedStart; j < alignedSize; j += RhsPacketSize)
499 _EIGEN_ACCUMULATE_PACKETS(Aligned, Unaligned, Aligned);
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)
540 _EIGEN_ACCUMULATE_PACKETS(Aligned, Unaligned, Unaligned);
541 break;
542 }
543 default:
544 for (Index j = alignedStart; j < alignedSize; j += RhsPacketSize)
545 _EIGEN_ACCUMULATE_PACKETS(Unaligned, Unaligned, Unaligned);
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 {
576 EIGEN_ALIGN_MAX ResScalar tmp0 = ResScalar(0);
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
Eigen::Solve
Pseudo expression representing a solving operation.
Definition Solve.h:63
Eigen::Unaligned
@ Unaligned
Data pointer has no specific alignment.
Definition Constants.h:228
Eigen::Aligned
@ Aligned
Definition Constants.h:235
Eigen::ColMajor
@ ColMajor
Storage order is column major (see TopicStorageOrders).
Definition Constants.h:320
Eigen::RowMajor
@ RowMajor
Storage order is row major (see TopicStorageOrders).
Definition Constants.h:322
Eigen::internal::general_matrix_vector_product
Definition BlasUtil.h:40
Eigen::internal::packet_traits
Definition GenericPacketMath.h:90
Eigen::internal::scalar_product_traits
Definition Meta.h:331
Eigen::internal::true_type
Definition Meta.h:30
Generated on Mon Sep 15 2025 12:11:54 for Medial Code Documentation by doxygen 1.9.8