25 #ifndef EIGEN_TRIANGULAR_MATRIX_MATRIX_H
26 #define EIGEN_TRIANGULAR_MATRIX_MATRIX_H
59 template <
typename Scalar,
typename Index,
60 int Mode,
bool LhsIsTriangular,
61 int LhsStorageOrder,
bool ConjugateLhs,
62 int RhsStorageOrder,
bool ConjugateRhs,
66 template <
typename Scalar,
typename Index,
67 int Mode,
bool LhsIsTriangular,
68 int LhsStorageOrder,
bool ConjugateLhs,
69 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
71 LhsStorageOrder,ConjugateLhs,
72 RhsStorageOrder,ConjugateRhs,
RowMajor,Version>
75 Index rows, Index cols, Index depth,
76 const Scalar* lhs, Index lhsStride,
77 const Scalar* rhs, Index rhsStride,
78 Scalar* res, Index resStride,
86 LhsStorageOrder==RowMajor ?
ColMajor : RowMajor,
89 ::run(cols, rows, depth, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha);
94 template <
typename Scalar,
typename Index,
int Mode,
95 int LhsStorageOrder,
bool ConjugateLhs,
96 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
98 LhsStorageOrder,ConjugateLhs,
99 RhsStorageOrder,ConjugateRhs,
ColMajor,Version>
102 typedef gebp_traits<Scalar,Scalar> Traits;
110 Index _rows, Index _cols, Index _depth,
111 const Scalar* _lhs, Index lhsStride,
112 const Scalar* _rhs, Index rhsStride,
113 Scalar* res, Index resStride,
117 Index diagSize = (std::min)(_rows,_depth);
118 Index rows = IsLower ? _rows : diagSize;
119 Index depth = IsLower ? diagSize : _depth;
122 const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
123 const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
128 computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
129 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
130 std::size_t sizeB = sizeW + kc*cols;
133 Scalar* blockB = allocatedBlockB + sizeW;
135 Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
136 triangularBuffer.setZero();
138 triangularBuffer.diagonal().setZero();
140 triangularBuffer.diagonal().setOnes();
142 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
143 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
144 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
146 for(Index k2=IsLower ? depth : 0;
147 IsLower ? k2>0 : k2<depth;
148 IsLower ? k2-=kc : k2+=kc)
150 Index actual_kc = (std::min)(IsLower ? k2 : depth-k2, kc);
151 Index actual_k2 = IsLower ? k2-actual_kc : k2;
154 if((!IsLower)&&(k2<rows)&&(k2+actual_kc>rows))
157 k2 = k2+actual_kc-kc;
160 pack_rhs(blockB, &rhs(actual_k2,0), rhsStride, actual_kc, cols);
168 if(IsLower || actual_k2<rows)
171 for (Index k1=0; k1<actual_kc; k1+=SmallPanelWidth)
173 Index actualPanelWidth = std::min<Index>(actual_kc-k1, SmallPanelWidth);
174 Index lengthTarget = IsLower ? actual_kc-k1-actualPanelWidth : k1;
175 Index startBlock = actual_k2+k1;
176 Index blockBOffset = k1;
181 for (Index k=0;k<actualPanelWidth;++k)
184 triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k);
185 for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i)
186 triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k);
188 pack_lhs(blockA, triangularBuffer.data(), triangularBuffer.outerStride(), actualPanelWidth, actualPanelWidth);
190 gebp_kernel(res+startBlock, resStride, blockA, blockB, actualPanelWidth, actualPanelWidth, cols, alpha,
191 actualPanelWidth, actual_kc, 0, blockBOffset);
196 Index startTarget = IsLower ? actual_k2+k1+actualPanelWidth : actual_k2;
198 pack_lhs(blockA, &lhs(startTarget,startBlock), lhsStride, actualPanelWidth, lengthTarget);
200 gebp_kernel(res+startTarget, resStride, blockA, blockB, lengthTarget, actualPanelWidth, cols, alpha,
201 actualPanelWidth, actual_kc, 0, blockBOffset);
207 Index start = IsLower ? k2 : 0;
208 Index end = IsLower ? rows : (std::min)(actual_k2,rows);
209 for(Index i2=start; i2<end; i2+=mc)
211 const Index actual_mc = (std::min)(i2+mc,end)-i2;
212 gemm_pack_lhs<Scalar, Index, Traits::mr,Traits::LhsProgress, LhsStorageOrder,false>()
213 (blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
215 gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
223 template <
typename Scalar,
typename Index,
int Mode,
224 int LhsStorageOrder,
bool ConjugateLhs,
225 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
227 LhsStorageOrder,ConjugateLhs,
228 RhsStorageOrder,ConjugateRhs,
ColMajor,Version>
230 typedef gebp_traits<Scalar,Scalar> Traits;
234 SetDiag = (Mode&(ZeroDiag|
UnitDiag)) ? 0 : 1
238 Index _rows, Index _cols, Index _depth,
239 const Scalar* _lhs, Index lhsStride,
240 const Scalar* _rhs, Index rhsStride,
241 Scalar* res, Index resStride,
245 Index diagSize = (std::min)(_cols,_depth);
247 Index depth = IsLower ? _depth : diagSize;
248 Index cols = IsLower ? diagSize : _cols;
250 const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
251 const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
256 computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
258 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
259 std::size_t sizeB = sizeW + kc*cols;
262 Scalar* blockB = allocatedBlockB + sizeW;
264 Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
265 triangularBuffer.setZero();
266 if((Mode&ZeroDiag)==ZeroDiag)
267 triangularBuffer.diagonal().setZero();
269 triangularBuffer.diagonal().setOnes();
271 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
272 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
273 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
274 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder,false,true> pack_rhs_panel;
276 for(Index k2=IsLower ? 0 : depth;
277 IsLower ? k2<depth : k2>0;
278 IsLower ? k2+=kc : k2-=kc)
280 Index actual_kc = (std::min)(IsLower ? depth-k2 : k2, kc);
281 Index actual_k2 = IsLower ? k2 : k2-actual_kc;
284 if(IsLower && (k2<cols) && (actual_k2+actual_kc>cols))
287 k2 = actual_k2 + actual_kc - kc;
291 Index rs = IsLower ? (std::min)(cols,actual_k2) : cols - k2;
293 Index ts = (IsLower && actual_k2>=cols) ? 0 : actual_kc;
295 Scalar* geb = blockB+ts*ts;
297 pack_rhs(geb, &rhs(actual_k2,IsLower ? 0 : k2), rhsStride, actual_kc, rs);
302 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
304 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
305 Index actual_j2 = actual_k2 + j2;
306 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
307 Index panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
309 pack_rhs_panel(blockB+j2*actual_kc,
310 &rhs(actual_k2+panelOffset, actual_j2), rhsStride,
311 panelLength, actualPanelWidth,
312 actual_kc, panelOffset);
315 for (Index j=0;j<actualPanelWidth;++j)
318 triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j);
319 for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k)
320 triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j);
323 pack_rhs_panel(blockB+j2*actual_kc,
324 triangularBuffer.data(), triangularBuffer.outerStride(),
325 actualPanelWidth, actualPanelWidth,
330 for (Index i2=0; i2<rows; i2+=mc)
332 const Index actual_mc = (std::min)(mc,rows-i2);
333 pack_lhs(blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
338 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
340 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
341 Index panelLength = IsLower ? actual_kc-j2 : j2+actualPanelWidth;
342 Index blockOffset = IsLower ? j2 : 0;
344 gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride,
345 blockA, blockB+j2*actual_kc,
346 actual_mc, panelLength, actualPanelWidth,
348 actual_kc, actual_kc,
349 blockOffset, blockOffset,
353 gebp_kernel(res+i2+(IsLower ? 0 : k2)*resStride, resStride,
354 blockA, geb, actual_mc, actual_kc, rs,
356 -1, -1, 0, 0, allocatedBlockB);
366 template<
int Mode,
bool LhsIsTriangular,
typename Lhs,
typename Rhs>
367 struct traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
368 :
traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>, Lhs, Rhs> >
373 template<
int Mode,
bool LhsIsTriangular,
typename Lhs,
typename Rhs>
374 struct TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
375 :
public ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>, Lhs, Rhs >
379 TriangularProduct(const Lhs& lhs, const Rhs& rhs) :
Base(lhs,rhs) {}
381 template<
typename Dest>
void scaleAndAddTo(Dest& dst,
Scalar alpha)
const
383 typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(m_lhs);
384 typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(m_rhs);
386 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
387 * RhsBlasTraits::extractScalarFactor(m_rhs);
389 internal::product_triangular_matrix_matrix<
Scalar,
Index,
390 Mode, LhsIsTriangular,
395 lhs.rows(), rhs.cols(), lhs.cols(),
396 &lhs.coeffRef(0,0), lhs.outerStride(),
397 &rhs.coeffRef(0,0), rhs.outerStride(),
398 &dst.coeffRef(0,0), dst.outerStride(),
406 #endif // EIGEN_TRIANGULAR_MATRIX_MATRIX_H