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