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-2016 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
15namespace internal {
16
17enum GEMVPacketSizeType {
18 GEMVPacketFull = 0,
19 GEMVPacketHalf,
20 GEMVPacketQuarter
21};
22
23template <int N, typename T1, typename T2, typename T3>
24struct gemv_packet_cond { typedef T3 type; };
25
26template <typename T1, typename T2, typename T3>
27struct gemv_packet_cond<GEMVPacketFull, T1, T2, T3> { typedef T1 type; };
28
29template <typename T1, typename T2, typename T3>
30struct gemv_packet_cond<GEMVPacketHalf, T1, T2, T3> { typedef T2 type; };
31
32template<typename LhsScalar, typename RhsScalar, int _PacketSize=GEMVPacketFull>
33class gemv_traits
34{
35 typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
36
37#define PACKET_DECL_COND_PREFIX(prefix, name, packet_size) \
38 typedef typename gemv_packet_cond<packet_size, \
39 typename packet_traits<name ## Scalar>::type, \
40 typename packet_traits<name ## Scalar>::half, \
41 typename unpacket_traits<typename packet_traits<name ## Scalar>::half>::half>::type \
42 prefix ## name ## Packet
43
44 PACKET_DECL_COND_PREFIX(_, Lhs, _PacketSize);
45 PACKET_DECL_COND_PREFIX(_, Rhs, _PacketSize);
46 PACKET_DECL_COND_PREFIX(_, Res, _PacketSize);
47#undef PACKET_DECL_COND_PREFIX
48
49public:
50 enum {
51 Vectorizable = unpacket_traits<_LhsPacket>::vectorizable &&
52 unpacket_traits<_RhsPacket>::vectorizable &&
53 int(unpacket_traits<_LhsPacket>::size)==int(unpacket_traits<_RhsPacket>::size),
54 LhsPacketSize = Vectorizable ? unpacket_traits<_LhsPacket>::size : 1,
55 RhsPacketSize = Vectorizable ? unpacket_traits<_RhsPacket>::size : 1,
56 ResPacketSize = Vectorizable ? unpacket_traits<_ResPacket>::size : 1
57 };
58
59 typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket;
60 typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket;
61 typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket;
62};
63
64
65/* Optimized col-major matrix * vector product:
66 * This algorithm processes the matrix per vertical panels,
67 * which are then processed horizontaly per chunck of 8*PacketSize x 1 vertical segments.
68 *
69 * Mixing type logic: C += alpha * A * B
70 * | A | B |alpha| comments
71 * |real |cplx |cplx | no vectorization
72 * |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization
73 * |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp
74 * |cplx |real |real | optimal case, vectorization possible via real-cplx mul
75 *
76 * The same reasoning apply for the transposed case.
77 */
78template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
79struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>
80{
81 typedef gemv_traits<LhsScalar,RhsScalar> Traits;
82 typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketHalf> HalfTraits;
83 typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketQuarter> QuarterTraits;
84
85 typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
86
87 typedef typename Traits::LhsPacket LhsPacket;
88 typedef typename Traits::RhsPacket RhsPacket;
89 typedef typename Traits::ResPacket ResPacket;
90
91 typedef typename HalfTraits::LhsPacket LhsPacketHalf;
92 typedef typename HalfTraits::RhsPacket RhsPacketHalf;
93 typedef typename HalfTraits::ResPacket ResPacketHalf;
94
95 typedef typename QuarterTraits::LhsPacket LhsPacketQuarter;
96 typedef typename QuarterTraits::RhsPacket RhsPacketQuarter;
97 typedef typename QuarterTraits::ResPacket ResPacketQuarter;
98
99EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run(
100 Index rows, Index cols,
101 const LhsMapper& lhs,
102 const RhsMapper& rhs,
103 ResScalar* res, Index resIncr,
104 RhsScalar alpha);
105};
106
107template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
108EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run(
109 Index rows, Index cols,
110 const LhsMapper& alhs,
111 const RhsMapper& rhs,
112 ResScalar* res, Index resIncr,
113 RhsScalar alpha)
114{
115 EIGEN_UNUSED_VARIABLE(resIncr);
116 eigen_internal_assert(resIncr==1);
117
118 // The following copy tells the compiler that lhs's attributes are not modified outside this function
119 // This helps GCC to generate propoer code.
120 LhsMapper lhs(alhs);
121
122 conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
123 conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj;
124 conj_helper<LhsPacketHalf,RhsPacketHalf,ConjugateLhs,ConjugateRhs> pcj_half;
125 conj_helper<LhsPacketQuarter,RhsPacketQuarter,ConjugateLhs,ConjugateRhs> pcj_quarter;
126
127 const Index lhsStride = lhs.stride();
128 // TODO: for padded aligned inputs, we could enable aligned reads
129 enum { LhsAlignment = Unaligned,
130 ResPacketSize = Traits::ResPacketSize,
131 ResPacketSizeHalf = HalfTraits::ResPacketSize,
132 ResPacketSizeQuarter = QuarterTraits::ResPacketSize,
133 LhsPacketSize = Traits::LhsPacketSize,
134 HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize,
135 HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf
136 };
137
138 const Index n8 = rows-8*ResPacketSize+1;
139 const Index n4 = rows-4*ResPacketSize+1;
140 const Index n3 = rows-3*ResPacketSize+1;
141 const Index n2 = rows-2*ResPacketSize+1;
142 const Index n1 = rows-1*ResPacketSize+1;
143 const Index n_half = rows-1*ResPacketSizeHalf+1;
144 const Index n_quarter = rows-1*ResPacketSizeQuarter+1;
145
146 // TODO: improve the following heuristic:
147 const Index block_cols = cols<128 ? cols : (lhsStride*sizeof(LhsScalar)<32000?16:4);
148 ResPacket palpha = pset1<ResPacket>(alpha);
149 ResPacketHalf palpha_half = pset1<ResPacketHalf>(alpha);
150 ResPacketQuarter palpha_quarter = pset1<ResPacketQuarter>(alpha);
151
152 for(Index j2=0; j2<cols; j2+=block_cols)
153 {
154 Index jend = numext::mini(j2+block_cols,cols);
155 Index i=0;
156 for(; i<n8; i+=ResPacketSize*8)
157 {
158 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
159 c1 = pset1<ResPacket>(ResScalar(0)),
160 c2 = pset1<ResPacket>(ResScalar(0)),
161 c3 = pset1<ResPacket>(ResScalar(0)),
162 c4 = pset1<ResPacket>(ResScalar(0)),
163 c5 = pset1<ResPacket>(ResScalar(0)),
164 c6 = pset1<ResPacket>(ResScalar(0)),
165 c7 = pset1<ResPacket>(ResScalar(0));
166
167 for(Index j=j2; j<jend; j+=1)
168 {
169 RhsPacket b0 = pset1<RhsPacket>(rhs(j,0));
170 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0);
171 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1);
172 c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2);
173 c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*3,j),b0,c3);
174 c4 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*4,j),b0,c4);
175 c5 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*5,j),b0,c5);
176 c6 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*6,j),b0,c6);
177 c7 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*7,j),b0,c7);
178 }
179 pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0)));
180 pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1)));
181 pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2)));
182 pstoreu(res+i+ResPacketSize*3, pmadd(c3,palpha,ploadu<ResPacket>(res+i+ResPacketSize*3)));
183 pstoreu(res+i+ResPacketSize*4, pmadd(c4,palpha,ploadu<ResPacket>(res+i+ResPacketSize*4)));
184 pstoreu(res+i+ResPacketSize*5, pmadd(c5,palpha,ploadu<ResPacket>(res+i+ResPacketSize*5)));
185 pstoreu(res+i+ResPacketSize*6, pmadd(c6,palpha,ploadu<ResPacket>(res+i+ResPacketSize*6)));
186 pstoreu(res+i+ResPacketSize*7, pmadd(c7,palpha,ploadu<ResPacket>(res+i+ResPacketSize*7)));
187 }
188 if(i<n4)
189 {
190 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
191 c1 = pset1<ResPacket>(ResScalar(0)),
192 c2 = pset1<ResPacket>(ResScalar(0)),
193 c3 = pset1<ResPacket>(ResScalar(0));
194
195 for(Index j=j2; j<jend; j+=1)
196 {
197 RhsPacket b0 = pset1<RhsPacket>(rhs(j,0));
198 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0);
199 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1);
200 c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2);
201 c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*3,j),b0,c3);
202 }
203 pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0)));
204 pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1)));
205 pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2)));
206 pstoreu(res+i+ResPacketSize*3, pmadd(c3,palpha,ploadu<ResPacket>(res+i+ResPacketSize*3)));
207
208 i+=ResPacketSize*4;
209 }
210 if(i<n3)
211 {
212 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
213 c1 = pset1<ResPacket>(ResScalar(0)),
214 c2 = pset1<ResPacket>(ResScalar(0));
215
216 for(Index j=j2; j<jend; j+=1)
217 {
218 RhsPacket b0 = pset1<RhsPacket>(rhs(j,0));
219 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0);
220 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1);
221 c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*2,j),b0,c2);
222 }
223 pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0)));
224 pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1)));
225 pstoreu(res+i+ResPacketSize*2, pmadd(c2,palpha,ploadu<ResPacket>(res+i+ResPacketSize*2)));
226
227 i+=ResPacketSize*3;
228 }
229 if(i<n2)
230 {
231 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
232 c1 = pset1<ResPacket>(ResScalar(0));
233
234 for(Index j=j2; j<jend; j+=1)
235 {
236 RhsPacket b0 = pset1<RhsPacket>(rhs(j,0));
237 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*0,j),b0,c0);
238 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+LhsPacketSize*1,j),b0,c1);
239 }
240 pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0)));
241 pstoreu(res+i+ResPacketSize*1, pmadd(c1,palpha,ploadu<ResPacket>(res+i+ResPacketSize*1)));
242 i+=ResPacketSize*2;
243 }
244 if(i<n1)
245 {
246 ResPacket c0 = pset1<ResPacket>(ResScalar(0));
247 for(Index j=j2; j<jend; j+=1)
248 {
249 RhsPacket b0 = pset1<RhsPacket>(rhs(j,0));
250 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0);
251 }
252 pstoreu(res+i+ResPacketSize*0, pmadd(c0,palpha,ploadu<ResPacket>(res+i+ResPacketSize*0)));
253 i+=ResPacketSize;
254 }
255 if(HasHalf && i<n_half)
256 {
257 ResPacketHalf c0 = pset1<ResPacketHalf>(ResScalar(0));
258 for(Index j=j2; j<jend; j+=1)
259 {
260 RhsPacketHalf b0 = pset1<RhsPacketHalf>(rhs(j,0));
261 c0 = pcj_half.pmadd(lhs.template load<LhsPacketHalf,LhsAlignment>(i+0,j),b0,c0);
262 }
263 pstoreu(res+i+ResPacketSizeHalf*0, pmadd(c0,palpha_half,ploadu<ResPacketHalf>(res+i+ResPacketSizeHalf*0)));
264 i+=ResPacketSizeHalf;
265 }
266 if(HasQuarter && i<n_quarter)
267 {
268 ResPacketQuarter c0 = pset1<ResPacketQuarter>(ResScalar(0));
269 for(Index j=j2; j<jend; j+=1)
270 {
271 RhsPacketQuarter b0 = pset1<RhsPacketQuarter>(rhs(j,0));
272 c0 = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter,LhsAlignment>(i+0,j),b0,c0);
273 }
274 pstoreu(res+i+ResPacketSizeQuarter*0, pmadd(c0,palpha_quarter,ploadu<ResPacketQuarter>(res+i+ResPacketSizeQuarter*0)));
275 i+=ResPacketSizeQuarter;
276 }
277 for(;i<rows;++i)
278 {
279 ResScalar c0(0);
280 for(Index j=j2; j<jend; j+=1)
281 c0 += cj.pmul(lhs(i,j), rhs(j,0));
282 res[i] += alpha*c0;
283 }
284 }
285}
286
287/* Optimized row-major matrix * vector product:
288 * This algorithm processes 4 rows at once that allows to both reduce
289 * the number of load/stores of the result by a factor 4 and to reduce
290 * the instruction dependency. Moreover, we know that all bands have the
291 * same alignment pattern.
292 *
293 * Mixing type logic:
294 * - alpha is always a complex (or converted to a complex)
295 * - no vectorization
296 */
297template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
298struct general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>
299{
300 typedef gemv_traits<LhsScalar,RhsScalar> Traits;
301 typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketHalf> HalfTraits;
302 typedef gemv_traits<LhsScalar,RhsScalar,GEMVPacketQuarter> QuarterTraits;
303
304 typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
305
306 typedef typename Traits::LhsPacket LhsPacket;
307 typedef typename Traits::RhsPacket RhsPacket;
308 typedef typename Traits::ResPacket ResPacket;
309
310 typedef typename HalfTraits::LhsPacket LhsPacketHalf;
311 typedef typename HalfTraits::RhsPacket RhsPacketHalf;
312 typedef typename HalfTraits::ResPacket ResPacketHalf;
313
314 typedef typename QuarterTraits::LhsPacket LhsPacketQuarter;
315 typedef typename QuarterTraits::RhsPacket RhsPacketQuarter;
316 typedef typename QuarterTraits::ResPacket ResPacketQuarter;
317
318EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE static void run(
319 Index rows, Index cols,
320 const LhsMapper& lhs,
321 const RhsMapper& rhs,
322 ResScalar* res, Index resIncr,
323 ResScalar alpha);
324};
325
326template<typename Index, typename LhsScalar, typename LhsMapper, bool ConjugateLhs, typename RhsScalar, typename RhsMapper, bool ConjugateRhs, int Version>
327EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,RowMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run(
328 Index rows, Index cols,
329 const LhsMapper& alhs,
330 const RhsMapper& rhs,
331 ResScalar* res, Index resIncr,
332 ResScalar alpha)
333{
334 // The following copy tells the compiler that lhs's attributes are not modified outside this function
335 // This helps GCC to generate propoer code.
336 LhsMapper lhs(alhs);
337
338 eigen_internal_assert(rhs.stride()==1);
339 conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj;
340 conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj;
341 conj_helper<LhsPacketHalf,RhsPacketHalf,ConjugateLhs,ConjugateRhs> pcj_half;
342 conj_helper<LhsPacketQuarter,RhsPacketQuarter,ConjugateLhs,ConjugateRhs> pcj_quarter;
343
344 // TODO: fine tune the following heuristic. The rationale is that if the matrix is very large,
345 // processing 8 rows at once might be counter productive wrt cache.
346 const Index n8 = lhs.stride()*sizeof(LhsScalar)>32000 ? 0 : rows-7;
347 const Index n4 = rows-3;
348 const Index n2 = rows-1;
349
350 // TODO: for padded aligned inputs, we could enable aligned reads
351 enum { LhsAlignment = Unaligned,
352 ResPacketSize = Traits::ResPacketSize,
353 ResPacketSizeHalf = HalfTraits::ResPacketSize,
354 ResPacketSizeQuarter = QuarterTraits::ResPacketSize,
355 LhsPacketSize = Traits::LhsPacketSize,
356 LhsPacketSizeHalf = HalfTraits::LhsPacketSize,
357 LhsPacketSizeQuarter = QuarterTraits::LhsPacketSize,
358 HasHalf = (int)ResPacketSizeHalf < (int)ResPacketSize,
359 HasQuarter = (int)ResPacketSizeQuarter < (int)ResPacketSizeHalf
360 };
361
362 typedef typename make_unsigned<Index>::type UnsignedIndex;
363 const Index fullColBlockEnd = LhsPacketSize * (UnsignedIndex(cols) / LhsPacketSize);
364 const Index halfColBlockEnd = LhsPacketSizeHalf * (UnsignedIndex(cols) / LhsPacketSizeHalf);
365 const Index quarterColBlockEnd = LhsPacketSizeQuarter * (UnsignedIndex(cols) / LhsPacketSizeQuarter);
366
367 Index i=0;
368 for(; i<n8; i+=8)
369 {
370 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
371 c1 = pset1<ResPacket>(ResScalar(0)),
372 c2 = pset1<ResPacket>(ResScalar(0)),
373 c3 = pset1<ResPacket>(ResScalar(0)),
374 c4 = pset1<ResPacket>(ResScalar(0)),
375 c5 = pset1<ResPacket>(ResScalar(0)),
376 c6 = pset1<ResPacket>(ResScalar(0)),
377 c7 = pset1<ResPacket>(ResScalar(0));
378
379 for (Index j = 0; j < fullColBlockEnd; j += LhsPacketSize)
380 {
381 RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0);
382
383 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0);
384 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1);
385 c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+2,j),b0,c2);
386 c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+3,j),b0,c3);
387 c4 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+4,j),b0,c4);
388 c5 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+5,j),b0,c5);
389 c6 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+6,j),b0,c6);
390 c7 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+7,j),b0,c7);
391 }
392 ResScalar cc0 = predux(c0);
393 ResScalar cc1 = predux(c1);
394 ResScalar cc2 = predux(c2);
395 ResScalar cc3 = predux(c3);
396 ResScalar cc4 = predux(c4);
397 ResScalar cc5 = predux(c5);
398 ResScalar cc6 = predux(c6);
399 ResScalar cc7 = predux(c7);
400
401 for (Index j = fullColBlockEnd; j < cols; ++j)
402 {
403 RhsScalar b0 = rhs(j,0);
404
405 cc0 += cj.pmul(lhs(i+0,j), b0);
406 cc1 += cj.pmul(lhs(i+1,j), b0);
407 cc2 += cj.pmul(lhs(i+2,j), b0);
408 cc3 += cj.pmul(lhs(i+3,j), b0);
409 cc4 += cj.pmul(lhs(i+4,j), b0);
410 cc5 += cj.pmul(lhs(i+5,j), b0);
411 cc6 += cj.pmul(lhs(i+6,j), b0);
412 cc7 += cj.pmul(lhs(i+7,j), b0);
413 }
414 res[(i+0)*resIncr] += alpha*cc0;
415 res[(i+1)*resIncr] += alpha*cc1;
416 res[(i+2)*resIncr] += alpha*cc2;
417 res[(i+3)*resIncr] += alpha*cc3;
418 res[(i+4)*resIncr] += alpha*cc4;
419 res[(i+5)*resIncr] += alpha*cc5;
420 res[(i+6)*resIncr] += alpha*cc6;
421 res[(i+7)*resIncr] += alpha*cc7;
422 }
423 for(; i<n4; i+=4)
424 {
425 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
426 c1 = pset1<ResPacket>(ResScalar(0)),
427 c2 = pset1<ResPacket>(ResScalar(0)),
428 c3 = pset1<ResPacket>(ResScalar(0));
429
430 for (Index j = 0; j < fullColBlockEnd; j += LhsPacketSize)
431 {
432 RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0);
433
434 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0);
435 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1);
436 c2 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+2,j),b0,c2);
437 c3 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+3,j),b0,c3);
438 }
439 ResScalar cc0 = predux(c0);
440 ResScalar cc1 = predux(c1);
441 ResScalar cc2 = predux(c2);
442 ResScalar cc3 = predux(c3);
443
444 for(Index j = fullColBlockEnd; j < cols; ++j)
445 {
446 RhsScalar b0 = rhs(j,0);
447
448 cc0 += cj.pmul(lhs(i+0,j), b0);
449 cc1 += cj.pmul(lhs(i+1,j), b0);
450 cc2 += cj.pmul(lhs(i+2,j), b0);
451 cc3 += cj.pmul(lhs(i+3,j), b0);
452 }
453 res[(i+0)*resIncr] += alpha*cc0;
454 res[(i+1)*resIncr] += alpha*cc1;
455 res[(i+2)*resIncr] += alpha*cc2;
456 res[(i+3)*resIncr] += alpha*cc3;
457 }
458 for(; i<n2; i+=2)
459 {
460 ResPacket c0 = pset1<ResPacket>(ResScalar(0)),
461 c1 = pset1<ResPacket>(ResScalar(0));
462
463 for (Index j = 0; j < fullColBlockEnd; j += LhsPacketSize)
464 {
465 RhsPacket b0 = rhs.template load<RhsPacket, Unaligned>(j,0);
466
467 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+0,j),b0,c0);
468 c1 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i+1,j),b0,c1);
469 }
470 ResScalar cc0 = predux(c0);
471 ResScalar cc1 = predux(c1);
472
473 for(Index j = fullColBlockEnd; j < cols; ++j)
474 {
475 RhsScalar b0 = rhs(j,0);
476
477 cc0 += cj.pmul(lhs(i+0,j), b0);
478 cc1 += cj.pmul(lhs(i+1,j), b0);
479 }
480 res[(i+0)*resIncr] += alpha*cc0;
481 res[(i+1)*resIncr] += alpha*cc1;
482 }
483 for(; i<rows; ++i)
484 {
485 ResPacket c0 = pset1<ResPacket>(ResScalar(0));
486 ResPacketHalf c0_h = pset1<ResPacketHalf>(ResScalar(0));
487 ResPacketQuarter c0_q = pset1<ResPacketQuarter>(ResScalar(0));
488
489 for (Index j = 0; j < fullColBlockEnd; j += LhsPacketSize)
490 {
491 RhsPacket b0 = rhs.template load<RhsPacket,Unaligned>(j,0);
492 c0 = pcj.pmadd(lhs.template load<LhsPacket,LhsAlignment>(i,j),b0,c0);
493 }
494 ResScalar cc0 = predux(c0);
495 if (HasHalf) {
496 for (Index j = fullColBlockEnd; j < halfColBlockEnd; j += LhsPacketSizeHalf)
497 {
498 RhsPacketHalf b0 = rhs.template load<RhsPacketHalf,Unaligned>(j,0);
499 c0_h = pcj_half.pmadd(lhs.template load<LhsPacketHalf,LhsAlignment>(i,j),b0,c0_h);
500 }
501 cc0 += predux(c0_h);
502 }
503 if (HasQuarter) {
504 for (Index j = halfColBlockEnd; j < quarterColBlockEnd; j += LhsPacketSizeQuarter)
505 {
506 RhsPacketQuarter b0 = rhs.template load<RhsPacketQuarter,Unaligned>(j,0);
507 c0_q = pcj_quarter.pmadd(lhs.template load<LhsPacketQuarter,LhsAlignment>(i,j),b0,c0_q);
508 }
509 cc0 += predux(c0_q);
510 }
511 for (Index j = quarterColBlockEnd; j < cols; ++j)
512 {
513 cc0 += cj.pmul(lhs(i,j), rhs(j,0));
514 }
515 res[i*resIncr] += alpha*cc0;
516 }
517}
518
519} // end namespace internal
520
521} // end namespace Eigen
522
523#endif // EIGEN_GENERAL_MATRIX_VECTOR_H
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::internal::gemv_traits
Definition GeneralMatrixVector.h:34
Eigen::Unaligned
@ Unaligned
Data pointer has no specific alignment.
Definition Constants.h:233
Eigen::ColMajor
@ ColMajor
Storage order is column major (see TopicStorageOrders).
Definition Constants.h:319
Eigen::RowMajor
@ RowMajor
Storage order is row major (see TopicStorageOrders).
Definition Constants.h:321
Eigen
Namespace containing all symbols from the Eigen library.
Definition LDLT.h:16
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition Meta.h:74
Eigen::ScalarBinaryOpTraits
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition XprHelper.h:806
Eigen::internal::gemv_packet_cond
Definition GeneralMatrixVector.h:24
Eigen::internal::general_matrix_vector_product
Definition BlasUtil.h:40
Eigen::internal::true_type
Definition Meta.h:96
Eigen::internal::unpacket_traits
Definition GenericPacketMath.h:133
Generated on Thu Dec 18 2025 18:56:20 for Medial Code Documentation by doxygen 1.9.8