10 #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_H 11 #define EIGEN_SELFADJOINT_MATRIX_MATRIX_H 18 template<
typename Scalar,
typename Index,
int Pack1,
int Pack2_dummy,
int StorageOrder>
21 template<
int BlockRows>
inline 27 blockA[count++] = lhs(i+
w,k);
30 for(
Index k=i; k<i+BlockRows; k++)
38 blockA[count++] = lhs(i+
w, k);
53 const Index peeled_mc3 = Pack1>=3*PacketSize ? (rows/(3*PacketSize))*(3*PacketSize) : 0;
54 const Index peeled_mc2 = Pack1>=2*PacketSize ? peeled_mc3+((rows-peeled_mc3)/(2*PacketSize))*(2*PacketSize) : 0;
55 const Index peeled_mc1 = Pack1>=1*PacketSize ? (rows/(1*PacketSize))*(1*PacketSize) : 0;
57 if(Pack1>=3*PacketSize)
58 for(
Index i=0;
i<peeled_mc3;
i+=3*PacketSize)
59 pack<3*PacketSize>(blockA, lhs, cols,
i, count);
61 if(Pack1>=2*PacketSize)
62 for(
Index i=peeled_mc3;
i<peeled_mc2;
i+=2*PacketSize)
63 pack<2*PacketSize>(blockA, lhs, cols,
i, count);
65 if(Pack1>=1*PacketSize)
66 for(
Index i=peeled_mc2;
i<peeled_mc1;
i+=1*PacketSize)
67 pack<1*PacketSize>(blockA, lhs, cols,
i, count);
73 blockA[count++] = lhs(i, k);
83 template<
typename Scalar,
typename Index,
int nr,
int StorageOrder>
92 Index packet_cols8 = nr>=8 ? (cols/8) * 8 : 0;
93 Index packet_cols4 = nr>=4 ? (cols/4) * 4 : 0;
96 for(
Index j2=0; j2<k2; j2+=nr)
98 for(
Index k=k2; k<end_k; k++)
100 blockB[count+0] = rhs(k,j2+0);
101 blockB[count+1] = rhs(k,j2+1);
104 blockB[count+2] = rhs(k,j2+2);
105 blockB[count+3] = rhs(k,j2+3);
109 blockB[count+4] = rhs(k,j2+4);
110 blockB[count+5] = rhs(k,j2+5);
111 blockB[count+6] = rhs(k,j2+6);
112 blockB[count+7] = rhs(k,j2+7);
122 for(
Index j2=k2; j2<end8; j2+=8)
126 for(
Index k=k2; k<j2; k++)
140 for(
Index k=j2; k<j2+8; k++)
144 blockB[count+
w] = rhs(k,j2+
w);
155 for(
Index k=j2+8; k<end_k; k++)
157 blockB[count+0] = rhs(k,j2+0);
158 blockB[count+1] = rhs(k,j2+1);
159 blockB[count+2] = rhs(k,j2+2);
160 blockB[count+3] = rhs(k,j2+3);
161 blockB[count+4] = rhs(k,j2+4);
162 blockB[count+5] = rhs(k,j2+5);
163 blockB[count+6] = rhs(k,j2+6);
164 blockB[count+7] = rhs(k,j2+7);
171 for(
Index j2=end8; j2<(
std::min)(k2+rows,packet_cols4); j2+=4)
175 for(
Index k=k2; k<j2; k++)
185 for(
Index k=j2; k<j2+4; k++)
189 blockB[count+
w] = rhs(k,j2+
w);
200 for(
Index k=j2+4; k<end_k; k++)
202 blockB[count+0] = rhs(k,j2+0);
203 blockB[count+1] = rhs(k,j2+1);
204 blockB[count+2] = rhs(k,j2+2);
205 blockB[count+3] = rhs(k,j2+3);
214 for(
Index j2=k2+rows; j2<packet_cols8; j2+=8)
216 for(
Index k=k2; k<end_k; k++)
232 for(
Index j2=(
std::max)(packet_cols8,k2+rows); j2<packet_cols4; j2+=4)
234 for(
Index k=k2; k<end_k; k++)
246 for(
Index j2=packet_cols4; j2<
cols; ++j2)
250 for(
Index k=k2; k<half; k++)
256 if(half==j2 && half<k2+rows)
267 blockB[count] = rhs(k,j2);
278 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
279 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs,
284 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
285 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs>
291 const Scalar* lhs, Index lhsStride,
292 const Scalar* rhs, Index rhsStride,
302 ::run(cols, rows, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha, blocking);
307 int LhsStorageOrder,
bool ConjugateLhs,
308 int RhsStorageOrder,
bool ConjugateRhs>
314 const Scalar* _lhs, Index lhsStride,
315 const Scalar* _rhs, Index rhsStride,
321 int LhsStorageOrder,
bool ConjugateLhs,
322 int RhsStorageOrder,
bool ConjugateRhs>
325 const Scalar* _lhs, Index lhsStride,
326 const Scalar* _rhs, Index rhsStride,
327 Scalar* _res, Index resStride,
338 LhsMapper lhs(_lhs,lhsStride);
339 LhsTransposeMapper lhs_transpose(_lhs,lhsStride);
340 RhsMapper rhs(_rhs,rhsStride);
341 ResMapper
res(_res, resStride);
343 Index kc = blocking.
kc();
357 for(Index k2=0; k2<
size; k2+=kc)
359 const Index actual_kc = (
std::min)(k2+kc,size)-k2;
364 pack_rhs(blockB, rhs.getSubMapper(k2,0), actual_kc,
cols);
370 for(Index i2=0; i2<k2; i2+=mc)
372 const Index actual_mc = (
std::min)(i2+mc,k2)-i2;
374 pack_lhs_transposed(blockA, lhs_transpose.getSubMapper(i2, k2), actual_kc, actual_mc);
376 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc,
cols,
alpha);
380 const Index actual_mc = (
std::min)(k2+kc,size)-k2;
382 pack_lhs(blockA, &lhs(k2,k2), lhsStride, actual_kc, actual_mc);
384 gebp_kernel(res.getSubMapper(k2, 0), blockA, blockB, actual_mc, actual_kc,
cols,
alpha);
387 for(Index i2=k2+kc; i2<
size; i2+=mc)
389 const Index actual_mc = (
std::min)(i2+mc,size)-i2;
391 (blockA, lhs.getSubMapper(i2, k2), actual_kc, actual_mc);
393 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc,
cols,
alpha);
400 int LhsStorageOrder,
bool ConjugateLhs,
401 int RhsStorageOrder,
bool ConjugateRhs>
407 const Scalar* _lhs, Index lhsStride,
408 const Scalar* _rhs, Index rhsStride,
409 Scalar*
res, Index resStride,
414 int LhsStorageOrder,
bool ConjugateLhs,
415 int RhsStorageOrder,
bool ConjugateRhs>
418 const Scalar* _lhs, Index lhsStride,
419 const Scalar* _rhs, Index rhsStride,
420 Scalar* _res, Index resStride,
429 LhsMapper lhs(_lhs,lhsStride);
430 ResMapper
res(_res,resStride);
432 Index kc = blocking.
kc();
443 for(Index k2=0; k2<
size; k2+=kc)
445 const Index actual_kc = (
std::min)(k2+kc,size)-k2;
447 pack_rhs(blockB, _rhs, rhsStride, actual_kc, cols, k2);
450 for(Index i2=0; i2<
rows; i2+=mc)
452 const Index actual_mc = (
std::min)(i2+mc,rows)-i2;
453 pack_lhs(blockA, lhs.getSubMapper(i2, k2), actual_kc, actual_mc);
455 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc,
cols,
alpha);
468 template<
typename Lhs,
int LhsMode,
typename Rhs,
int RhsMode>
485 template<
typename Dest>
486 static void run(Dest &dst,
const Lhs &a_lhs,
const Rhs &a_rhs,
const Scalar&
alpha)
488 eigen_assert(dst.rows()==a_lhs.rows() && dst.cols()==a_rhs.cols());
493 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(a_lhs)
494 * RhsBlasTraits::extractScalarFactor(a_rhs);
497 Lhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxColsAtCompileTime,1> BlockingType;
499 BlockingType blocking(lhs.rows(), rhs.cols(), lhs.cols(), 1,
false);
508 lhs.rows(), rhs.cols(),
509 &lhs.coeffRef(0,0), lhs.outerStride(),
510 &rhs.coeffRef(0,0), rhs.outerStride(),
511 &dst.coeffRef(0,0), dst.outerStride(),
512 actualAlpha, blocking
521 #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_H LhsBlasTraits::DirectLinearAccessType ActualLhsType
void operator()(Scalar *blockB, const Scalar *_rhs, Index rhsStride, Index rows, Index cols, Index k2)
#define EIGEN_STRONG_INLINE
Expression of the product of two arbitrary matrices or vectors.
Product< Lhs, Rhs >::Scalar Scalar
Namespace containing all symbols from the Eigen library.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
#define EIGEN_LOGICAL_XOR(a, b)
internal::blas_traits< Lhs > LhsBlasTraits
const unsigned int RowMajorBit
internal::blas_traits< Rhs > RhsBlasTraits
#define EIGEN_DONT_INLINE
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
void pack(Scalar *blockA, const const_blas_data_mapper< Scalar, Index, StorageOrder > &lhs, Index cols, Index i, Index &count)
RhsBlasTraits::DirectLinearAccessType ActualRhsType
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
void operator()(Scalar *blockA, const Scalar *_lhs, Index lhsStride, Index cols, Index rows)
void run(Expr &expr, Dev &dev)
static void run(Dest &dst, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar &alpha)
static EIGEN_STRONG_INLINE void run(Index rows, Index cols, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsStride, Scalar *res, Index resStride, const Scalar &alpha, level3_blocking< Scalar, Scalar > &blocking)
ScalarWithExceptions conj(const ScalarWithExceptions &x)