10 #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_H
11 #define EIGEN_SELFADJOINT_MATRIX_MATRIX_H
18 template<
typename Scalar,
typename Index,
int Pack1,
int Pack2,
int StorageOrder>
21 template<
int BlockRows>
inline
22 void pack(Scalar* blockA,
const const_blas_data_mapper<Scalar,Index,StorageOrder>& lhs, Index cols, Index i, Index& count)
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++)
32 for(Index w=0; w<h; w++)
33 blockA[count++] =
conj(lhs(k, i+w));
35 blockA[count++] =
real(lhs(k,k));
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++] =
conj(lhs(k, i+w));
46 void operator()(Scalar* blockA,
const Scalar* _lhs, Index lhsStride, Index cols, Index rows)
48 const_blas_data_mapper<Scalar,Index,StorageOrder> lhs(_lhs,lhsStride);
50 Index peeled_mc = (rows/Pack1)*Pack1;
51 for(Index i=0; i<peeled_mc; i+=Pack1)
53 pack<Pack1>(blockA, lhs, cols, i, count);
56 if(rows-peeled_mc>=Pack2)
58 pack<Pack2>(blockA, lhs, cols, peeled_mc, count);
63 for(Index i=peeled_mc; i<rows; i++)
65 for(Index k=0; k<i; k++)
66 blockA[count++] = lhs(i, k);
68 blockA[count++] =
real(lhs(i, i));
70 for(Index k=i+1; k<cols; k++)
71 blockA[count++] =
conj(lhs(k, i));
76 template<
typename Scalar,
typename Index,
int nr,
int StorageOrder>
79 enum { PacketSize = packet_traits<Scalar>::size };
80 void operator()(Scalar* blockB,
const Scalar* _rhs, Index rhsStride, Index rows, Index cols, Index k2)
82 Index end_k = k2 + rows;
84 const_blas_data_mapper<Scalar,Index,StorageOrder> rhs(_rhs,rhsStride);
85 Index packet_cols = (cols/nr)*nr;
88 for(Index j2=0; j2<k2; j2+=nr)
90 for(Index k=k2; k<end_k; k++)
92 blockB[count+0] = rhs(k,j2+0);
93 blockB[count+1] = rhs(k,j2+1);
96 blockB[count+2] = rhs(k,j2+2);
97 blockB[count+3] = rhs(k,j2+3);
104 for(Index j2=k2; j2<(std::min)(k2+rows,packet_cols); j2+=nr)
108 for(Index k=k2; k<j2; k++)
110 blockB[count+0] =
conj(rhs(j2+0,k));
111 blockB[count+1] =
conj(rhs(j2+1,k));
114 blockB[count+2] =
conj(rhs(j2+2,k));
115 blockB[count+3] =
conj(rhs(j2+3,k));
121 for(Index k=j2; k<j2+nr; k++)
124 for (Index w=0 ; w<h; ++w)
125 blockB[count+w] = rhs(k,j2+w);
127 blockB[count+h] =
real(rhs(k,k));
130 for (Index w=h+1 ; w<nr; ++w)
131 blockB[count+w] =
conj(rhs(j2+w,k));
136 for(Index k=j2+nr; k<end_k; k++)
138 blockB[count+0] = rhs(k,j2+0);
139 blockB[count+1] = rhs(k,j2+1);
142 blockB[count+2] = rhs(k,j2+2);
143 blockB[count+3] = rhs(k,j2+3);
150 for(Index j2=k2+rows; j2<packet_cols; j2+=nr)
152 for(Index k=k2; k<end_k; k++)
154 blockB[count+0] =
conj(rhs(j2+0,k));
155 blockB[count+1] =
conj(rhs(j2+1,k));
158 blockB[count+2] =
conj(rhs(j2+2,k));
159 blockB[count+3] =
conj(rhs(j2+3,k));
166 for(Index j2=packet_cols; j2<cols; ++j2)
169 Index half = (std::min)(end_k,j2);
170 for(Index k=k2; k<half; k++)
172 blockB[count] =
conj(rhs(j2,k));
176 if(half==j2 && half<k2+rows)
178 blockB[count] =
real(rhs(j2,j2));
185 for(Index k=half+1; k<k2+rows; k++)
187 blockB[count] = rhs(k,j2);
197 template <
typename Scalar,
typename Index,
198 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
199 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs,
201 struct product_selfadjoint_matrix;
203 template <
typename Scalar,
typename Index,
204 int LhsStorageOrder,
bool LhsSelfAdjoint,
bool ConjugateLhs,
205 int RhsStorageOrder,
bool RhsSelfAdjoint,
bool ConjugateRhs>
206 struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,LhsSelfAdjoint,ConjugateLhs, RhsStorageOrder,RhsSelfAdjoint,ConjugateRhs,
RowMajor>
209 static EIGEN_STRONG_INLINE
void run(
210 Index rows, Index cols,
211 const Scalar* lhs, Index lhsStride,
212 const Scalar* rhs, Index rhsStride,
213 Scalar* res, Index resStride,
216 product_selfadjoint_matrix<Scalar, Index,
218 RhsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsSelfAdjoint,ConjugateRhs),
220 LhsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsSelfAdjoint,ConjugateLhs),
222 ::run(cols, rows, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha);
226 template <
typename Scalar,
typename Index,
227 int LhsStorageOrder,
bool ConjugateLhs,
228 int RhsStorageOrder,
bool ConjugateRhs>
229 struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs, RhsStorageOrder,false,ConjugateRhs,
ColMajor>
232 static EIGEN_DONT_INLINE
void run(
233 Index rows, Index cols,
234 const Scalar* _lhs, Index lhsStride,
235 const Scalar* _rhs, Index rhsStride,
236 Scalar* res, Index resStride,
241 const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
242 const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
244 typedef gebp_traits<Scalar,Scalar> Traits;
249 computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
251 kc = (std::min)(kc,mc);
253 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
254 std::size_t sizeB = sizeW + kc*cols;
255 ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
256 ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
257 Scalar* blockB = allocatedBlockB + sizeW;
259 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
260 symm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
261 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
262 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder==RowMajor?ColMajor:RowMajor, true> pack_lhs_transposed;
264 for(Index k2=0; k2<size; k2+=kc)
266 const Index actual_kc = (std::min)(k2+kc,size)-k2;
271 pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, cols);
277 for(Index i2=0; i2<k2; i2+=mc)
279 const Index actual_mc = (std::min)(i2+mc,k2)-i2;
281 pack_lhs_transposed(blockA, &lhs(k2, i2), lhsStride, actual_kc, actual_mc);
283 gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
287 const Index actual_mc = (std::min)(k2+kc,size)-k2;
289 pack_lhs(blockA, &lhs(k2,k2), lhsStride, actual_kc, actual_mc);
291 gebp_kernel(res+k2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
294 for(Index i2=k2+kc; i2<size; i2+=mc)
296 const Index actual_mc = (std::min)(i2+mc,size)-i2;
297 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder,false>()
298 (blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
300 gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
307 template <
typename Scalar,
typename Index,
308 int LhsStorageOrder,
bool ConjugateLhs,
309 int RhsStorageOrder,
bool ConjugateRhs>
310 struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLhs, RhsStorageOrder,true,ConjugateRhs,
ColMajor>
313 static EIGEN_DONT_INLINE
void run(
314 Index rows, Index cols,
315 const Scalar* _lhs, Index lhsStride,
316 const Scalar* _rhs, Index rhsStride,
317 Scalar* res, Index resStride,
322 const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
324 typedef gebp_traits<Scalar,Scalar> Traits;
329 computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
330 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
331 std::size_t sizeB = sizeW + kc*cols;
332 ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
333 ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
334 Scalar* blockB = allocatedBlockB + sizeW;
336 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
337 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
338 symm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
340 for(Index k2=0; k2<size; k2+=kc)
342 const Index actual_kc = (std::min)(k2+kc,size)-k2;
344 pack_rhs(blockB, _rhs, rhsStride, actual_kc, cols, k2);
347 for(Index i2=0; i2<rows; i2+=mc)
349 const Index actual_mc = (std::min)(i2+mc,rows)-i2;
350 pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
352 gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
365 template<
typename Lhs,
int LhsMode,
typename Rhs,
int RhsMode>
366 struct traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
367 : traits<ProductBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>, Lhs, Rhs> >
371 template<
typename Lhs,
int LhsMode,
typename Rhs,
int RhsMode>
372 struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
373 :
public ProductBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>, Lhs, Rhs >
375 EIGEN_PRODUCT_PUBLIC_INTERFACE(SelfadjointProductMatrix)
377 SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
383 RhsIsSelfAdjoint = (RhsMode&
SelfAdjoint)==SelfAdjoint
386 template<
typename Dest>
void scaleAndAddTo(Dest& dst, Scalar alpha)
const
388 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
390 typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(m_lhs);
391 typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(m_rhs);
393 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
394 * RhsBlasTraits::extractScalarFactor(m_rhs);
396 internal::product_selfadjoint_matrix<Scalar, Index,
397 EIGEN_LOGICAL_XOR(LhsIsUpper,
399 NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsIsUpper,
bool(LhsBlasTraits::NeedToConjugate)),
400 EIGEN_LOGICAL_XOR(RhsIsUpper,
402 NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsIsUpper,
bool(RhsBlasTraits::NeedToConjugate)),
405 lhs.rows(), rhs.cols(),
406 &lhs.coeffRef(0,0), lhs.outerStride(),
407 &rhs.coeffRef(0,0), rhs.outerStride(),
408 &dst.coeffRef(0,0), dst.outerStride(),
416 #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_H