Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
GeneralMatrixVector.h
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_GENERAL_MATRIX_VECTOR_H 00011 #define EIGEN_GENERAL_MATRIX_VECTOR_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00017 /* Optimized col-major matrix * vector product: 00018 * This algorithm processes 4 columns at onces that allows to both reduce 00019 * the number of load/stores of the result by a factor 4 and to reduce 00020 * the instruction dependency. Moreover, we know that all bands have the 00021 * same alignment pattern. 00022 * 00023 * Mixing type logic: C += alpha * A * B 00024 * | A | B |alpha| comments 00025 * |real |cplx |cplx | no vectorization 00026 * |real |cplx |real | alpha is converted to a cplx when calling the run function, no vectorization 00027 * |cplx |real |cplx | invalid, the caller has to do tmp: = A * B; C += alpha*tmp 00028 * |cplx |real |real | optimal case, vectorization possible via real-cplx mul 00029 */ 00030 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version> 00031 struct general_matrix_vector_product<Index,LhsScalar,ColMajor,ConjugateLhs,RhsScalar,ConjugateRhs,Version> 00032 { 00033 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; 00034 00035 enum { 00036 Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable 00037 && int(packet_traits<LhsScalar>::size)==int(packet_traits<RhsScalar>::size), 00038 LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1, 00039 RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1, 00040 ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1 00041 }; 00042 00043 typedef typename packet_traits<LhsScalar>::type _LhsPacket; 00044 typedef typename packet_traits<RhsScalar>::type _RhsPacket; 00045 typedef typename packet_traits<ResScalar>::type _ResPacket; 00046 00047 typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket; 00048 typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket; 00049 typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket; 00050 00051 EIGEN_DONT_INLINE static void run( 00052 Index rows, Index cols, 00053 const LhsScalar* lhs, Index lhsStride, 00054 const RhsScalar* rhs, Index rhsIncr, 00055 ResScalar* res, Index resIncr, RhsScalar alpha); 00056 }; 00057 00058 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version> 00059 EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,ColMajor,ConjugateLhs,RhsScalar,ConjugateRhs,Version>::run( 00060 Index rows, Index cols, 00061 const LhsScalar* lhs, Index lhsStride, 00062 const RhsScalar* rhs, Index rhsIncr, 00063 ResScalar* res, Index resIncr, RhsScalar alpha) 00064 { 00065 EIGEN_UNUSED_VARIABLE(resIncr) 00066 eigen_internal_assert(resIncr==1); 00067 #ifdef _EIGEN_ACCUMULATE_PACKETS 00068 #error _EIGEN_ACCUMULATE_PACKETS has already been defined 00069 #endif 00070 #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) \ 00071 pstore(&res[j], \ 00072 padd(pload<ResPacket>(&res[j]), \ 00073 padd( \ 00074 padd(pcj.pmul(EIGEN_CAT(ploa , A0)<LhsPacket>(&lhs0[j]), ptmp0), \ 00075 pcj.pmul(EIGEN_CAT(ploa , A13)<LhsPacket>(&lhs1[j]), ptmp1)), \ 00076 padd(pcj.pmul(EIGEN_CAT(ploa , A2)<LhsPacket>(&lhs2[j]), ptmp2), \ 00077 pcj.pmul(EIGEN_CAT(ploa , A13)<LhsPacket>(&lhs3[j]), ptmp3)) ))) 00078 00079 conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj; 00080 conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj; 00081 if(ConjugateRhs) 00082 alpha = numext::conj(alpha); 00083 00084 enum { AllAligned = 0, EvenAligned, FirstAligned, NoneAligned }; 00085 const Index columnsAtOnce = 4; 00086 const Index peels = 2; 00087 const Index LhsPacketAlignedMask = LhsPacketSize-1; 00088 const Index ResPacketAlignedMask = ResPacketSize-1; 00089 // const Index PeelAlignedMask = ResPacketSize*peels-1; 00090 const Index size = rows; 00091 00092 // How many coeffs of the result do we have to skip to be aligned. 00093 // Here we assume data are at least aligned on the base scalar type. 00094 Index alignedStart = internal::first_aligned(res,size); 00095 Index alignedSize = ResPacketSize>1 ? alignedStart + ((size-alignedStart) & ~ResPacketAlignedMask) : 0; 00096 const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1; 00097 00098 const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0; 00099 Index alignmentPattern = alignmentStep==0 ? AllAligned 00100 : alignmentStep==(LhsPacketSize/2) ? EvenAligned 00101 : FirstAligned; 00102 00103 // we cannot assume the first element is aligned because of sub-matrices 00104 const Index lhsAlignmentOffset = internal::first_aligned(lhs,size); 00105 00106 // find how many columns do we have to skip to be aligned with the result (if possible) 00107 Index skipColumns = 0; 00108 // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats) 00109 if( (size_t(lhs)%sizeof(LhsScalar)) || (size_t(res)%sizeof(ResScalar)) ) 00110 { 00111 alignedSize = 0; 00112 alignedStart = 0; 00113 } 00114 else if (LhsPacketSize>1) 00115 { 00116 eigen_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0 || size<LhsPacketSize); 00117 00118 while (skipColumns<LhsPacketSize && 00119 alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%LhsPacketSize)) 00120 ++skipColumns; 00121 if (skipColumns==LhsPacketSize) 00122 { 00123 // nothing can be aligned, no need to skip any column 00124 alignmentPattern = NoneAligned; 00125 skipColumns = 0; 00126 } 00127 else 00128 { 00129 skipColumns = (std::min)(skipColumns,cols); 00130 // note that the skiped columns are processed later. 00131 } 00132 00133 eigen_internal_assert( (alignmentPattern==NoneAligned) 00134 || (skipColumns + columnsAtOnce >= cols) 00135 || LhsPacketSize > size 00136 || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(LhsPacket))==0); 00137 } 00138 else if(Vectorizable) 00139 { 00140 alignedStart = 0; 00141 alignedSize = size; 00142 alignmentPattern = AllAligned; 00143 } 00144 00145 Index offset1 = (FirstAligned && alignmentStep==1?3:1); 00146 Index offset3 = (FirstAligned && alignmentStep==1?1:3); 00147 00148 Index columnBound = ((cols-skipColumns)/columnsAtOnce)*columnsAtOnce + skipColumns; 00149 for (Index i=skipColumns; i<columnBound; i+=columnsAtOnce) 00150 { 00151 RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs[i*rhsIncr]), 00152 ptmp1 = pset1<RhsPacket>(alpha*rhs[(i+offset1)*rhsIncr]), 00153 ptmp2 = pset1<RhsPacket>(alpha*rhs[(i+2)*rhsIncr]), 00154 ptmp3 = pset1<RhsPacket>(alpha*rhs[(i+offset3)*rhsIncr]); 00155 00156 // this helps a lot generating better binary code 00157 const LhsScalar *lhs0 = lhs + i*lhsStride, *lhs1 = lhs + (i+offset1)*lhsStride, 00158 *lhs2 = lhs + (i+2)*lhsStride, *lhs3 = lhs + (i+offset3)*lhsStride; 00159 00160 if (Vectorizable) 00161 { 00162 /* explicit vectorization */ 00163 // process initial unaligned coeffs 00164 for (Index j=0; j<alignedStart; ++j) 00165 { 00166 res[j] = cj.pmadd(lhs0[j], pfirst(ptmp0), res[j]); 00167 res[j] = cj.pmadd(lhs1[j], pfirst(ptmp1), res[j]); 00168 res[j] = cj.pmadd(lhs2[j], pfirst(ptmp2), res[j]); 00169 res[j] = cj.pmadd(lhs3[j], pfirst(ptmp3), res[j]); 00170 } 00171 00172 if (alignedSize>alignedStart) 00173 { 00174 switch(alignmentPattern) 00175 { 00176 case AllAligned: 00177 for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize) 00178 _EIGEN_ACCUMULATE_PACKETS(d,d,d); 00179 break; 00180 case EvenAligned: 00181 for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize) 00182 _EIGEN_ACCUMULATE_PACKETS(d,du,d); 00183 break; 00184 case FirstAligned: 00185 { 00186 Index j = alignedStart; 00187 if(peels>1) 00188 { 00189 LhsPacket A00, A01, A02, A03, A10, A11, A12, A13; 00190 ResPacket T0, T1; 00191 00192 A01 = pload<LhsPacket>(&lhs1[alignedStart-1]); 00193 A02 = pload<LhsPacket>(&lhs2[alignedStart-2]); 00194 A03 = pload<LhsPacket>(&lhs3[alignedStart-3]); 00195 00196 for (; j<peeledSize; j+=peels*ResPacketSize) 00197 { 00198 A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11); 00199 A12 = pload<LhsPacket>(&lhs2[j-2+LhsPacketSize]); palign<2>(A02,A12); 00200 A13 = pload<LhsPacket>(&lhs3[j-3+LhsPacketSize]); palign<3>(A03,A13); 00201 00202 A00 = pload<LhsPacket>(&lhs0[j]); 00203 A10 = pload<LhsPacket>(&lhs0[j+LhsPacketSize]); 00204 T0 = pcj.pmadd(A00, ptmp0, pload<ResPacket>(&res[j])); 00205 T1 = pcj.pmadd(A10, ptmp0, pload<ResPacket>(&res[j+ResPacketSize])); 00206 00207 T0 = pcj.pmadd(A01, ptmp1, T0); 00208 A01 = pload<LhsPacket>(&lhs1[j-1+2*LhsPacketSize]); palign<1>(A11,A01); 00209 T0 = pcj.pmadd(A02, ptmp2, T0); 00210 A02 = pload<LhsPacket>(&lhs2[j-2+2*LhsPacketSize]); palign<2>(A12,A02); 00211 T0 = pcj.pmadd(A03, ptmp3, T0); 00212 pstore(&res[j],T0); 00213 A03 = pload<LhsPacket>(&lhs3[j-3+2*LhsPacketSize]); palign<3>(A13,A03); 00214 T1 = pcj.pmadd(A11, ptmp1, T1); 00215 T1 = pcj.pmadd(A12, ptmp2, T1); 00216 T1 = pcj.pmadd(A13, ptmp3, T1); 00217 pstore(&res[j+ResPacketSize],T1); 00218 } 00219 } 00220 for (; j<alignedSize; j+=ResPacketSize) 00221 _EIGEN_ACCUMULATE_PACKETS(d,du,du); 00222 break; 00223 } 00224 default: 00225 for (Index j = alignedStart; j<alignedSize; j+=ResPacketSize) 00226 _EIGEN_ACCUMULATE_PACKETS(du,du,du); 00227 break; 00228 } 00229 } 00230 } // end explicit vectorization 00231 00232 /* process remaining coeffs (or all if there is no explicit vectorization) */ 00233 for (Index j=alignedSize; j<size; ++j) 00234 { 00235 res[j] = cj.pmadd(lhs0[j], pfirst(ptmp0), res[j]); 00236 res[j] = cj.pmadd(lhs1[j], pfirst(ptmp1), res[j]); 00237 res[j] = cj.pmadd(lhs2[j], pfirst(ptmp2), res[j]); 00238 res[j] = cj.pmadd(lhs3[j], pfirst(ptmp3), res[j]); 00239 } 00240 } 00241 00242 // process remaining first and last columns (at most columnsAtOnce-1) 00243 Index end = cols; 00244 Index start = columnBound; 00245 do 00246 { 00247 for (Index k=start; k<end; ++k) 00248 { 00249 RhsPacket ptmp0 = pset1<RhsPacket>(alpha*rhs[k*rhsIncr]); 00250 const LhsScalar* lhs0 = lhs + k*lhsStride; 00251 00252 if (Vectorizable) 00253 { 00254 /* explicit vectorization */ 00255 // process first unaligned result's coeffs 00256 for (Index j=0; j<alignedStart; ++j) 00257 res[j] += cj.pmul(lhs0[j], pfirst(ptmp0)); 00258 // process aligned result's coeffs 00259 if ((size_t(lhs0+alignedStart)%sizeof(LhsPacket))==0) 00260 for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize) 00261 pstore(&res[i], pcj.pmadd(pload<LhsPacket>(&lhs0[i]), ptmp0, pload<ResPacket>(&res[i]))); 00262 else 00263 for (Index i = alignedStart;i<alignedSize;i+=ResPacketSize) 00264 pstore(&res[i], pcj.pmadd(ploadu<LhsPacket>(&lhs0[i]), ptmp0, pload<ResPacket>(&res[i]))); 00265 } 00266 00267 // process remaining scalars (or all if no explicit vectorization) 00268 for (Index i=alignedSize; i<size; ++i) 00269 res[i] += cj.pmul(lhs0[i], pfirst(ptmp0)); 00270 } 00271 if (skipColumns) 00272 { 00273 start = 0; 00274 end = skipColumns; 00275 skipColumns = 0; 00276 } 00277 else 00278 break; 00279 } while(Vectorizable); 00280 #undef _EIGEN_ACCUMULATE_PACKETS 00281 } 00282 00283 /* Optimized row-major matrix * vector product: 00284 * This algorithm processes 4 rows at onces that allows to both reduce 00285 * the number of load/stores of the result by a factor 4 and to reduce 00286 * the instruction dependency. Moreover, we know that all bands have the 00287 * same alignment pattern. 00288 * 00289 * Mixing type logic: 00290 * - alpha is always a complex (or converted to a complex) 00291 * - no vectorization 00292 */ 00293 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version> 00294 struct general_matrix_vector_product<Index,LhsScalar,RowMajor,ConjugateLhs,RhsScalar,ConjugateRhs,Version> 00295 { 00296 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar; 00297 00298 enum { 00299 Vectorizable = packet_traits<LhsScalar>::Vectorizable && packet_traits<RhsScalar>::Vectorizable 00300 && int(packet_traits<LhsScalar>::size)==int(packet_traits<RhsScalar>::size), 00301 LhsPacketSize = Vectorizable ? packet_traits<LhsScalar>::size : 1, 00302 RhsPacketSize = Vectorizable ? packet_traits<RhsScalar>::size : 1, 00303 ResPacketSize = Vectorizable ? packet_traits<ResScalar>::size : 1 00304 }; 00305 00306 typedef typename packet_traits<LhsScalar>::type _LhsPacket; 00307 typedef typename packet_traits<RhsScalar>::type _RhsPacket; 00308 typedef typename packet_traits<ResScalar>::type _ResPacket; 00309 00310 typedef typename conditional<Vectorizable,_LhsPacket,LhsScalar>::type LhsPacket; 00311 typedef typename conditional<Vectorizable,_RhsPacket,RhsScalar>::type RhsPacket; 00312 typedef typename conditional<Vectorizable,_ResPacket,ResScalar>::type ResPacket; 00313 00314 EIGEN_DONT_INLINE static void run( 00315 Index rows, Index cols, 00316 const LhsScalar* lhs, Index lhsStride, 00317 const RhsScalar* rhs, Index rhsIncr, 00318 ResScalar* res, Index resIncr, 00319 ResScalar alpha); 00320 }; 00321 00322 template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version> 00323 EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,RowMajor,ConjugateLhs,RhsScalar,ConjugateRhs,Version>::run( 00324 Index rows, Index cols, 00325 const LhsScalar* lhs, Index lhsStride, 00326 const RhsScalar* rhs, Index rhsIncr, 00327 ResScalar* res, Index resIncr, 00328 ResScalar alpha) 00329 { 00330 EIGEN_UNUSED_VARIABLE(rhsIncr); 00331 eigen_internal_assert(rhsIncr==1); 00332 #ifdef _EIGEN_ACCUMULATE_PACKETS 00333 #error _EIGEN_ACCUMULATE_PACKETS has already been defined 00334 #endif 00335 00336 #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) {\ 00337 RhsPacket b = pload<RhsPacket>(&rhs[j]); \ 00338 ptmp0 = pcj.pmadd(EIGEN_CAT(ploa,A0) <LhsPacket>(&lhs0[j]), b, ptmp0); \ 00339 ptmp1 = pcj.pmadd(EIGEN_CAT(ploa,A13)<LhsPacket>(&lhs1[j]), b, ptmp1); \ 00340 ptmp2 = pcj.pmadd(EIGEN_CAT(ploa,A2) <LhsPacket>(&lhs2[j]), b, ptmp2); \ 00341 ptmp3 = pcj.pmadd(EIGEN_CAT(ploa,A13)<LhsPacket>(&lhs3[j]), b, ptmp3); } 00342 00343 conj_helper<LhsScalar,RhsScalar,ConjugateLhs,ConjugateRhs> cj; 00344 conj_helper<LhsPacket,RhsPacket,ConjugateLhs,ConjugateRhs> pcj; 00345 00346 enum { AllAligned=0, EvenAligned=1, FirstAligned=2, NoneAligned=3 }; 00347 const Index rowsAtOnce = 4; 00348 const Index peels = 2; 00349 const Index RhsPacketAlignedMask = RhsPacketSize-1; 00350 const Index LhsPacketAlignedMask = LhsPacketSize-1; 00351 // const Index PeelAlignedMask = RhsPacketSize*peels-1; 00352 const Index depth = cols; 00353 00354 // How many coeffs of the result do we have to skip to be aligned. 00355 // Here we assume data are at least aligned on the base scalar type 00356 // if that's not the case then vectorization is discarded, see below. 00357 Index alignedStart = internal::first_aligned(rhs, depth); 00358 Index alignedSize = RhsPacketSize>1 ? alignedStart + ((depth-alignedStart) & ~RhsPacketAlignedMask) : 0; 00359 const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1; 00360 00361 const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0; 00362 Index alignmentPattern = alignmentStep==0 ? AllAligned 00363 : alignmentStep==(LhsPacketSize/2) ? EvenAligned 00364 : FirstAligned; 00365 00366 // we cannot assume the first element is aligned because of sub-matrices 00367 const Index lhsAlignmentOffset = internal::first_aligned(lhs,depth); 00368 00369 // find how many rows do we have to skip to be aligned with rhs (if possible) 00370 Index skipRows = 0; 00371 // if the data cannot be aligned (TODO add some compile time tests when possible, e.g. for floats) 00372 if( (sizeof(LhsScalar)!=sizeof(RhsScalar)) || (size_t(lhs)%sizeof(LhsScalar)) || (size_t(rhs)%sizeof(RhsScalar)) ) 00373 { 00374 alignedSize = 0; 00375 alignedStart = 0; 00376 } 00377 else if (LhsPacketSize>1) 00378 { 00379 eigen_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(LhsPacket)==0 || depth<LhsPacketSize); 00380 00381 while (skipRows<LhsPacketSize && 00382 alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%LhsPacketSize)) 00383 ++skipRows; 00384 if (skipRows==LhsPacketSize) 00385 { 00386 // nothing can be aligned, no need to skip any column 00387 alignmentPattern = NoneAligned; 00388 skipRows = 0; 00389 } 00390 else 00391 { 00392 skipRows = (std::min)(skipRows,Index(rows)); 00393 // note that the skiped columns are processed later. 00394 } 00395 eigen_internal_assert( alignmentPattern==NoneAligned 00396 || LhsPacketSize==1 00397 || (skipRows + rowsAtOnce >= rows) 00398 || LhsPacketSize > depth 00399 || (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(LhsPacket))==0); 00400 } 00401 else if(Vectorizable) 00402 { 00403 alignedStart = 0; 00404 alignedSize = depth; 00405 alignmentPattern = AllAligned; 00406 } 00407 00408 Index offset1 = (FirstAligned && alignmentStep==1?3:1); 00409 Index offset3 = (FirstAligned && alignmentStep==1?1:3); 00410 00411 Index rowBound = ((rows-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows; 00412 for (Index i=skipRows; i<rowBound; i+=rowsAtOnce) 00413 { 00414 EIGEN_ALIGN16 ResScalar tmp0 = ResScalar(0); 00415 ResScalar tmp1 = ResScalar(0), tmp2 = ResScalar(0), tmp3 = ResScalar(0); 00416 00417 // this helps the compiler generating good binary code 00418 const LhsScalar *lhs0 = lhs + i*lhsStride, *lhs1 = lhs + (i+offset1)*lhsStride, 00419 *lhs2 = lhs + (i+2)*lhsStride, *lhs3 = lhs + (i+offset3)*lhsStride; 00420 00421 if (Vectorizable) 00422 { 00423 /* explicit vectorization */ 00424 ResPacket ptmp0 = pset1<ResPacket>(ResScalar(0)), ptmp1 = pset1<ResPacket>(ResScalar(0)), 00425 ptmp2 = pset1<ResPacket>(ResScalar(0)), ptmp3 = pset1<ResPacket>(ResScalar(0)); 00426 00427 // process initial unaligned coeffs 00428 // FIXME this loop get vectorized by the compiler ! 00429 for (Index j=0; j<alignedStart; ++j) 00430 { 00431 RhsScalar b = rhs[j]; 00432 tmp0 += cj.pmul(lhs0[j],b); tmp1 += cj.pmul(lhs1[j],b); 00433 tmp2 += cj.pmul(lhs2[j],b); tmp3 += cj.pmul(lhs3[j],b); 00434 } 00435 00436 if (alignedSize>alignedStart) 00437 { 00438 switch(alignmentPattern) 00439 { 00440 case AllAligned: 00441 for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize) 00442 _EIGEN_ACCUMULATE_PACKETS(d,d,d); 00443 break; 00444 case EvenAligned: 00445 for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize) 00446 _EIGEN_ACCUMULATE_PACKETS(d,du,d); 00447 break; 00448 case FirstAligned: 00449 { 00450 Index j = alignedStart; 00451 if (peels>1) 00452 { 00453 /* Here we proccess 4 rows with with two peeled iterations to hide 00454 * the overhead of unaligned loads. Moreover unaligned loads are handled 00455 * using special shift/move operations between the two aligned packets 00456 * overlaping the desired unaligned packet. This is *much* more efficient 00457 * than basic unaligned loads. 00458 */ 00459 LhsPacket A01, A02, A03, A11, A12, A13; 00460 A01 = pload<LhsPacket>(&lhs1[alignedStart-1]); 00461 A02 = pload<LhsPacket>(&lhs2[alignedStart-2]); 00462 A03 = pload<LhsPacket>(&lhs3[alignedStart-3]); 00463 00464 for (; j<peeledSize; j+=peels*RhsPacketSize) 00465 { 00466 RhsPacket b = pload<RhsPacket>(&rhs[j]); 00467 A11 = pload<LhsPacket>(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11); 00468 A12 = pload<LhsPacket>(&lhs2[j-2+LhsPacketSize]); palign<2>(A02,A12); 00469 A13 = pload<LhsPacket>(&lhs3[j-3+LhsPacketSize]); palign<3>(A03,A13); 00470 00471 ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j]), b, ptmp0); 00472 ptmp1 = pcj.pmadd(A01, b, ptmp1); 00473 A01 = pload<LhsPacket>(&lhs1[j-1+2*LhsPacketSize]); palign<1>(A11,A01); 00474 ptmp2 = pcj.pmadd(A02, b, ptmp2); 00475 A02 = pload<LhsPacket>(&lhs2[j-2+2*LhsPacketSize]); palign<2>(A12,A02); 00476 ptmp3 = pcj.pmadd(A03, b, ptmp3); 00477 A03 = pload<LhsPacket>(&lhs3[j-3+2*LhsPacketSize]); palign<3>(A13,A03); 00478 00479 b = pload<RhsPacket>(&rhs[j+RhsPacketSize]); 00480 ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j+LhsPacketSize]), b, ptmp0); 00481 ptmp1 = pcj.pmadd(A11, b, ptmp1); 00482 ptmp2 = pcj.pmadd(A12, b, ptmp2); 00483 ptmp3 = pcj.pmadd(A13, b, ptmp3); 00484 } 00485 } 00486 for (; j<alignedSize; j+=RhsPacketSize) 00487 _EIGEN_ACCUMULATE_PACKETS(d,du,du); 00488 break; 00489 } 00490 default: 00491 for (Index j = alignedStart; j<alignedSize; j+=RhsPacketSize) 00492 _EIGEN_ACCUMULATE_PACKETS(du,du,du); 00493 break; 00494 } 00495 tmp0 += predux(ptmp0); 00496 tmp1 += predux(ptmp1); 00497 tmp2 += predux(ptmp2); 00498 tmp3 += predux(ptmp3); 00499 } 00500 } // end explicit vectorization 00501 00502 // process remaining coeffs (or all if no explicit vectorization) 00503 // FIXME this loop get vectorized by the compiler ! 00504 for (Index j=alignedSize; j<depth; ++j) 00505 { 00506 RhsScalar b = rhs[j]; 00507 tmp0 += cj.pmul(lhs0[j],b); tmp1 += cj.pmul(lhs1[j],b); 00508 tmp2 += cj.pmul(lhs2[j],b); tmp3 += cj.pmul(lhs3[j],b); 00509 } 00510 res[i*resIncr] += alpha*tmp0; 00511 res[(i+offset1)*resIncr] += alpha*tmp1; 00512 res[(i+2)*resIncr] += alpha*tmp2; 00513 res[(i+offset3)*resIncr] += alpha*tmp3; 00514 } 00515 00516 // process remaining first and last rows (at most columnsAtOnce-1) 00517 Index end = rows; 00518 Index start = rowBound; 00519 do 00520 { 00521 for (Index i=start; i<end; ++i) 00522 { 00523 EIGEN_ALIGN16 ResScalar tmp0 = ResScalar(0); 00524 ResPacket ptmp0 = pset1<ResPacket>(tmp0); 00525 const LhsScalar* lhs0 = lhs + i*lhsStride; 00526 // process first unaligned result's coeffs 00527 // FIXME this loop get vectorized by the compiler ! 00528 for (Index j=0; j<alignedStart; ++j) 00529 tmp0 += cj.pmul(lhs0[j], rhs[j]); 00530 00531 if (alignedSize>alignedStart) 00532 { 00533 // process aligned rhs coeffs 00534 if ((size_t(lhs0+alignedStart)%sizeof(LhsPacket))==0) 00535 for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize) 00536 ptmp0 = pcj.pmadd(pload<LhsPacket>(&lhs0[j]), pload<RhsPacket>(&rhs[j]), ptmp0); 00537 else 00538 for (Index j = alignedStart;j<alignedSize;j+=RhsPacketSize) 00539 ptmp0 = pcj.pmadd(ploadu<LhsPacket>(&lhs0[j]), pload<RhsPacket>(&rhs[j]), ptmp0); 00540 tmp0 += predux(ptmp0); 00541 } 00542 00543 // process remaining scalars 00544 // FIXME this loop get vectorized by the compiler ! 00545 for (Index j=alignedSize; j<depth; ++j) 00546 tmp0 += cj.pmul(lhs0[j], rhs[j]); 00547 res[i*resIncr] += alpha*tmp0; 00548 } 00549 if (skipRows) 00550 { 00551 start = 0; 00552 end = skipRows; 00553 skipRows = 0; 00554 } 00555 else 00556 break; 00557 } while(Vectorizable); 00558 00559 #undef _EIGEN_ACCUMULATE_PACKETS 00560 } 00561 00562 } // end namespace internal 00563 00564 } // end namespace Eigen 00565 00566 #endif // EIGEN_GENERAL_MATRIX_VECTOR_H
Generated on Tue Jul 12 2022 17:46:54 by 1.7.2