10#ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_H
11#define EIGEN_SELFADJOINT_MATRIX_MATRIX_H
18template<
typename Scalar,
typename Index,
int Pack1,
int Pack2_dummy,
int StorageOrder>
21 template<
int BlockRows>
inline
25 for(
Index k=0; k<i; k++)
26 for(
Index w=0; w<BlockRows; w++)
27 blockA[
count++] = lhs(i+w,k);
30 for(
Index k=i; k<i+BlockRows; k++)
33 blockA[
count++] = numext::conj(lhs(k, i+w));
37 for(
Index w=
h+1; w<BlockRows; w++)
38 blockA[
count++] = lhs(i+w, k);
42 for(
Index k=i+BlockRows; k<cols; k++)
43 for(
Index w=0; w<BlockRows; w++)
44 blockA[
count++] = numext::conj(lhs(k, i+w));
53 HasHalf = (int)HalfPacketSize < (
int)PacketSize,
54 HasQuarter = (int)QuarterPacketSize < (
int)HalfPacketSize};
60 const Index peeled_mc3 = Pack1>=3*PacketSize ? (rows/(3*PacketSize))*(3*PacketSize) : 0;
61 const Index peeled_mc2 = Pack1>=2*PacketSize ? peeled_mc3+((rows-peeled_mc3)/(2*PacketSize))*(2*PacketSize) : 0;
62 const Index peeled_mc1 = Pack1>=1*PacketSize ? peeled_mc2+((rows-peeled_mc2)/(1*PacketSize))*(1*PacketSize) : 0;
63 const Index peeled_mc_half = Pack1>=HalfPacketSize ? peeled_mc1+((rows-peeled_mc1)/(HalfPacketSize))*(HalfPacketSize) : 0;
64 const Index peeled_mc_quarter = Pack1>=QuarterPacketSize ? peeled_mc_half+((rows-peeled_mc_half)/(QuarterPacketSize))*(QuarterPacketSize) : 0;
66 if(Pack1>=3*PacketSize)
67 for(
Index i=0; i<peeled_mc3; i+=3*PacketSize)
68 pack<3*PacketSize>(blockA, lhs, cols, i,
count);
70 if(Pack1>=2*PacketSize)
71 for(
Index i=peeled_mc3; i<peeled_mc2; i+=2*PacketSize)
72 pack<2*PacketSize>(blockA, lhs, cols, i,
count);
74 if(Pack1>=1*PacketSize)
75 for(
Index i=peeled_mc2; i<peeled_mc1; i+=1*PacketSize)
76 pack<1*PacketSize>(blockA, lhs, cols, i,
count);
78 if(HasHalf && Pack1>=HalfPacketSize)
79 for(
Index i=peeled_mc1; i<peeled_mc_half; i+=HalfPacketSize)
80 pack<HalfPacketSize>(blockA, lhs, cols, i,
count);
82 if(HasQuarter && Pack1>=QuarterPacketSize)
83 for(
Index i=peeled_mc_half; i<peeled_mc_quarter; i+=QuarterPacketSize)
84 pack<QuarterPacketSize>(blockA, lhs, cols, i,
count);
87 for(
Index i=peeled_mc_quarter; i<rows; i++)
89 for(
Index k=0; k<i; k++)
90 blockA[
count++] = lhs(i, k);
94 for(
Index k=i+1; k<cols; k++)
95 blockA[
count++] = numext::conj(lhs(k, i));
100template<
typename Scalar,
typename Index,
int nr,
int StorageOrder>
109 Index packet_cols8 = nr>=8 ? (cols/8) * 8 : 0;
110 Index packet_cols4 = nr>=4 ? (cols/4) * 4 : 0;
117 blockB[
count+0] = rhs(k,j2+0);
118 blockB[
count+1] = rhs(k,j2+1);
121 blockB[
count+2] = rhs(k,j2+2);
122 blockB[
count+3] = rhs(k,j2+3);
126 blockB[
count+4] = rhs(k,j2+4);
127 blockB[
count+5] = rhs(k,j2+5);
128 blockB[
count+6] = rhs(k,j2+6);
129 blockB[
count+7] = rhs(k,j2+7);
139 for(
Index j2=
k2; j2<end8; j2+=8)
145 blockB[
count+0] = numext::conj(rhs(j2+0,k));
146 blockB[
count+1] = numext::conj(rhs(j2+1,k));
147 blockB[
count+2] = numext::conj(rhs(j2+2,k));
148 blockB[
count+3] = numext::conj(rhs(j2+3,k));
149 blockB[
count+4] = numext::conj(rhs(j2+4,k));
150 blockB[
count+5] = numext::conj(rhs(j2+5,k));
151 blockB[
count+6] = numext::conj(rhs(j2+6,k));
152 blockB[
count+7] = numext::conj(rhs(j2+7,k));
157 for(
Index k=j2; k<j2+8; k++)
160 for (
Index w=0 ; w<
h; ++w)
161 blockB[
count+w] = rhs(k,j2+w);
166 for (
Index w=
h+1 ; w<8; ++w)
167 blockB[
count+w] = numext::conj(rhs(j2+w,k));
172 for(
Index k=j2+8; k<end_k; k++)
174 blockB[
count+0] = rhs(k,j2+0);
175 blockB[
count+1] = rhs(k,j2+1);
176 blockB[
count+2] = rhs(k,j2+2);
177 blockB[
count+3] = rhs(k,j2+3);
178 blockB[
count+4] = rhs(k,j2+4);
179 blockB[
count+5] = rhs(k,j2+5);
180 blockB[
count+6] = rhs(k,j2+6);
181 blockB[
count+7] = rhs(k,j2+7);
194 blockB[
count+0] = numext::conj(rhs(j2+0,k));
195 blockB[
count+1] = numext::conj(rhs(j2+1,k));
196 blockB[
count+2] = numext::conj(rhs(j2+2,k));
197 blockB[
count+3] = numext::conj(rhs(j2+3,k));
202 for(
Index k=j2; k<j2+4; k++)
205 for (
Index w=0 ; w<
h; ++w)
206 blockB[
count+w] = rhs(k,j2+w);
211 for (
Index w=
h+1 ; w<4; ++w)
212 blockB[
count+w] = numext::conj(rhs(j2+w,k));
217 for(
Index k=j2+4; k<end_k; k++)
219 blockB[
count+0] = rhs(k,j2+0);
220 blockB[
count+1] = rhs(k,j2+1);
221 blockB[
count+2] = rhs(k,j2+2);
222 blockB[
count+3] = rhs(k,j2+3);
231 for(
Index j2=
k2+rows; j2<packet_cols8; j2+=8)
235 blockB[
count+0] = numext::conj(rhs(j2+0,k));
236 blockB[
count+1] = numext::conj(rhs(j2+1,k));
237 blockB[
count+2] = numext::conj(rhs(j2+2,k));
238 blockB[
count+3] = numext::conj(rhs(j2+3,k));
239 blockB[
count+4] = numext::conj(rhs(j2+4,k));
240 blockB[
count+5] = numext::conj(rhs(j2+5,k));
241 blockB[
count+6] = numext::conj(rhs(j2+6,k));
242 blockB[
count+7] = numext::conj(rhs(j2+7,k));
253 blockB[
count+0] = numext::conj(rhs(j2+0,k));
254 blockB[
count+1] = numext::conj(rhs(j2+1,k));
255 blockB[
count+2] = numext::conj(rhs(j2+2,k));
256 blockB[
count+3] = numext::conj(rhs(j2+3,k));
263 for(
Index j2=packet_cols4; j2<cols; ++j2)
269 blockB[
count] = numext::conj(rhs(j2,k));
284 blockB[
count] = rhs(k,j2);
294template <
typename Scalar,
typename Index,
295 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
296 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs,
297 int ResStorageOrder,
int ResInnerStride>
300template <
typename Scalar,
typename Index,
301 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
302 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs,
309 const Scalar* lhs,
Index lhsStride,
310 const Scalar* rhs,
Index rhsStride,
320 ::run(cols, rows, rhs, rhsStride, lhs, lhsStride, res, resIncr, resStride, alpha, blocking);
324template <
typename Scalar,
typename Index,
325 int LhsStorageOrder,
bool ConjugateLhs,
326 int RhsStorageOrder,
bool ConjugateRhs,
333 const Scalar* _lhs,
Index lhsStride,
334 const Scalar* _rhs,
Index rhsStride,
339template <
typename Scalar,
typename Index,
340 int LhsStorageOrder,
bool ConjugateLhs,
341 int RhsStorageOrder,
bool ConjugateRhs,
345 const Scalar* _lhs,
Index lhsStride,
346 const Scalar* _rhs,
Index rhsStride,
358 LhsMapper lhs(_lhs,lhsStride);
359 LhsTransposeMapper lhs_transpose(_lhs,lhsStride);
360 RhsMapper rhs(_rhs,rhsStride);
361 ResMapper res(_res, resStride, resIncr);
367 std::size_t sizeA = kc*mc;
368 std::size_t sizeB = kc*cols;
384 pack_rhs(blockB, rhs.getSubMapper(
k2,0), actual_kc, cols);
394 pack_lhs_transposed(blockA, lhs_transpose.getSubMapper(i2,
k2), actual_kc, actual_mc);
396 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc, cols, alpha);
402 pack_lhs(blockA, &lhs(
k2,
k2), lhsStride, actual_kc, actual_mc);
404 gebp_kernel(res.getSubMapper(
k2, 0), blockA, blockB, actual_mc, actual_kc, cols, alpha);
411 (blockA, lhs.getSubMapper(i2,
k2), actual_kc, actual_mc);
413 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc, cols, alpha);
419template <
typename Scalar,
typename Index,
420 int LhsStorageOrder,
bool ConjugateLhs,
421 int RhsStorageOrder,
bool ConjugateRhs,
428 const Scalar* _lhs,
Index lhsStride,
429 const Scalar* _rhs,
Index rhsStride,
434template <
typename Scalar,
typename Index,
435 int LhsStorageOrder,
bool ConjugateLhs,
436 int RhsStorageOrder,
bool ConjugateRhs,
440 const Scalar* _lhs,
Index lhsStride,
441 const Scalar* _rhs,
Index rhsStride,
451 LhsMapper lhs(_lhs,lhsStride);
452 ResMapper res(_res,resStride, resIncr);
456 std::size_t sizeA = kc*mc;
457 std::size_t sizeB = kc*cols;
469 pack_rhs(blockB, _rhs, rhsStride, actual_kc, cols,
k2);
472 for(
Index i2=0; i2<rows; i2+=mc)
475 pack_lhs(blockA, lhs.getSubMapper(i2,
k2), actual_kc, actual_mc);
477 gebp_kernel(res.getSubMapper(i2, 0), blockA, blockB, actual_mc, actual_kc, cols, alpha);
490template<
typename Lhs,
int LhsMode,
typename Rhs,
int RhsMode>
507 template<
typename Dest>
508 static void run(Dest &dst,
const Lhs &a_lhs,
const Rhs &a_rhs,
const Scalar& alpha)
510 eigen_assert(dst.rows()==a_lhs.rows() && dst.cols()==a_rhs.cols());
515 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(a_lhs)
516 * RhsBlasTraits::extractScalarFactor(a_rhs);
519 Lhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxColsAtCompileTime,1> BlockingType;
521 BlockingType blocking(lhs.rows(), rhs.cols(), lhs.cols(), 1,
false);
529 Dest::InnerStrideAtCompileTime>
531 lhs.rows(), rhs.cols(),
532 &lhs.coeffRef(0,0), lhs.outerStride(),
533 &rhs.coeffRef(0,0), rhs.outerStride(),
534 &dst.coeffRef(0,0), dst.innerStride(), dst.outerStride(),
535 actualAlpha, blocking
EIGEN_DEVICE_FUNC RealReturnType real() const
Definition: CommonCwiseUnaryOps.h:100
#define EIGEN_LOGICAL_XOR(a, b)
Definition: Macros.h:1323
#define EIGEN_DONT_INLINE
Definition: Macros.h:950
#define eigen_assert(x)
Definition: Macros.h:1047
#define EIGEN_STRONG_INLINE
Definition: Macros.h:927
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:768
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:75
Definition: BlasUtil.h:270
Definition: BlasUtil.h:389
Definition: GeneralBlockPanelKernel.h:419
Definition: GeneralMatrixMatrix.h:248
Definition: GeneralMatrixMatrix.h:252
RhsScalar * blockB()
Definition: GeneralMatrixMatrix.h:275
Index kc() const
Definition: GeneralMatrixMatrix.h:272
Index mc() const
Definition: GeneralMatrixMatrix.h:270
LhsScalar * blockA()
Definition: GeneralMatrixMatrix.h:274
constexpr auto count() -> size_t
Definition: core.h:1204
@ SelfAdjoint
Used in BandMatrix and SelfAdjointView to indicate that the matrix is self-adjoint.
Definition: Constants.h:225
@ Lower
View matrix as a lower triangular matrix.
Definition: Constants.h:209
@ Upper
View matrix as an upper triangular matrix.
Definition: Constants.h:211
@ ColMajor
Storage order is column major (see TopicStorageOrders).
Definition: Constants.h:319
@ RowMajor
Storage order is row major (see TopicStorageOrders).
Definition: Constants.h:321
const unsigned int RowMajorBit
for a matrix, this means that the storage order is row-major.
Definition: Constants.h:66
constexpr common_t< T1, T2 > max(const T1 x, const T2 y) noexcept
Compile-time pairwise maximum function.
Definition: max.hpp:35
constexpr common_t< T1, T2 > min(const T1 x, const T2 y) noexcept
Compile-time pairwise minimum function.
Definition: min.hpp:35
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
Namespace containing all symbols from the Eigen library.
Definition: MatrixExponential.h:16
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Eigen_Colamd.h:50
static constexpr const unit_t< compound_unit< energy::joule, time::seconds > > h(6.626070040e-34)
Planck constant.
static constexpr uint64_t k2
Definition: Hashing.h:172
Holds information about the various numeric (i.e.
Definition: NumTraits.h:233
const T type
Definition: Meta.h:214
Definition: BlasUtil.h:403
Definition: GeneralBlockPanelKernel.h:1058
Definition: BlasUtil.h:28
Definition: BlasUtil.h:25
Definition: GenericPacketMath.h:107
static EIGEN_STRONG_INLINE void run(Index rows, Index cols, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsStride, Scalar *res, Index resIncr, Index resStride, const Scalar &alpha, level3_blocking< Scalar, Scalar > &blocking)
Definition: SelfadjointMatrixMatrix.h:307
Definition: SelfadjointMatrixMatrix.h:298
internal::blas_traits< Lhs > LhsBlasTraits
Definition: SelfadjointMatrixMatrix.h:495
Product< Lhs, Rhs >::Scalar Scalar
Definition: SelfadjointMatrixMatrix.h:493
RhsBlasTraits::DirectLinearAccessType ActualRhsType
Definition: SelfadjointMatrixMatrix.h:498
LhsBlasTraits::DirectLinearAccessType ActualLhsType
Definition: SelfadjointMatrixMatrix.h:496
static void run(Dest &dst, const Lhs &a_lhs, const Rhs &a_rhs, const Scalar &alpha)
Definition: SelfadjointMatrixMatrix.h:508
internal::blas_traits< Rhs > RhsBlasTraits
Definition: SelfadjointMatrixMatrix.h:497
Definition: ProductEvaluators.h:793
Definition: SelfadjointMatrixMatrix.h:20
void pack(Scalar *blockA, const const_blas_data_mapper< Scalar, Index, StorageOrder > &lhs, Index cols, Index i, Index &count)
Definition: SelfadjointMatrixMatrix.h:22
void operator()(Scalar *blockA, const Scalar *_lhs, Index lhsStride, Index cols, Index rows)
Definition: SelfadjointMatrixMatrix.h:46
Definition: SelfadjointMatrixMatrix.h:102
void operator()(Scalar *blockB, const Scalar *_rhs, Index rhsStride, Index rows, Index cols, Index k2)
Definition: SelfadjointMatrixMatrix.h:104
@ PacketSize
Definition: SelfadjointMatrixMatrix.h:103
Definition: ForwardDeclarations.h:17
Definition: GenericPacketMath.h:133