10 #ifndef EIGEN_TRIANGULAR_MATRIX_MATRIX_H 11 #define EIGEN_TRIANGULAR_MATRIX_MATRIX_H 44 template <
typename Scalar,
typename Index,
45 int Mode,
bool LhsIsTriangular,
46 int LhsStorageOrder,
bool ConjugateLhs,
47 int RhsStorageOrder,
bool ConjugateRhs,
51 template <
typename Scalar,
typename Index,
52 int Mode,
bool LhsIsTriangular,
53 int LhsStorageOrder,
bool ConjugateLhs,
54 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
56 LhsStorageOrder,ConjugateLhs,
57 RhsStorageOrder,ConjugateRhs,
RowMajor,Version>
60 Index rows, Index cols, Index depth,
61 const Scalar* lhs, Index lhsStride,
62 const Scalar*
rhs, Index rhsStride,
63 Scalar* res, Index resStride,
71 LhsStorageOrder==RowMajor ?
ColMajor : RowMajor,
74 ::run(cols, rows, depth, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha, blocking);
79 template <
typename Scalar,
typename Index,
int Mode,
80 int LhsStorageOrder,
bool ConjugateLhs,
81 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
83 LhsStorageOrder,ConjugateLhs,
84 RhsStorageOrder,ConjugateRhs,
ColMajor,Version>
95 Index _rows, Index _cols, Index _depth,
96 const Scalar* _lhs, Index lhsStride,
97 const Scalar* _rhs, Index rhsStride,
98 Scalar* res, Index resStride,
102 template <
typename Scalar,
typename Index,
int Mode,
103 int LhsStorageOrder,
bool ConjugateLhs,
104 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
106 LhsStorageOrder,ConjugateLhs,
107 RhsStorageOrder,ConjugateRhs,
ColMajor,Version>::run(
108 Index _rows, Index _cols, Index _depth,
109 const Scalar* _lhs, Index lhsStride,
110 const Scalar* _rhs, Index rhsStride,
111 Scalar* res, Index resStride,
115 Index diagSize = (std::min)(_rows,_depth);
116 Index rows = IsLower ? _rows : diagSize;
117 Index depth = IsLower ? diagSize : _depth;
123 Index kc = blocking.
kc();
124 Index mc = (std::min)(rows,blocking.
mc());
126 std::size_t sizeA = kc*mc;
127 std::size_t sizeB = kc*cols;
128 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
137 triangularBuffer.diagonal().setZero();
139 triangularBuffer.diagonal().setOnes();
145 for(Index k2=IsLower ? depth : 0;
146 IsLower ? k2>0 : k2<depth;
147 IsLower ? k2-=kc : k2+=kc)
149 Index actual_kc = (std::min)(IsLower ? k2 : depth-k2, kc);
150 Index actual_k2 = IsLower ? k2-actual_kc : k2;
153 if((!IsLower)&&(k2<rows)&&(k2+actual_kc>rows))
156 k2 = k2+actual_kc-kc;
159 pack_rhs(blockB, &
rhs(actual_k2,0), rhsStride, actual_kc, cols);
167 if(IsLower || actual_k2<rows)
170 for (Index k1=0; k1<actual_kc; k1+=SmallPanelWidth)
172 Index actualPanelWidth = std::min<Index>(actual_kc-k1, SmallPanelWidth);
173 Index lengthTarget = IsLower ? actual_kc-k1-actualPanelWidth : k1;
174 Index startBlock = actual_k2+k1;
175 Index blockBOffset = k1;
180 for (Index k=0;k<actualPanelWidth;++k)
183 triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k);
184 for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i)
185 triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k);
187 pack_lhs(blockA, triangularBuffer.data(), triangularBuffer.outerStride(), actualPanelWidth, actualPanelWidth);
189 gebp_kernel(res+startBlock, resStride, blockA, blockB, actualPanelWidth, actualPanelWidth, cols, alpha,
190 actualPanelWidth, actual_kc, 0, blockBOffset, blockW);
195 Index startTarget = IsLower ? actual_k2+k1+actualPanelWidth : actual_k2;
197 pack_lhs(blockA, &lhs(startTarget,startBlock), lhsStride, actualPanelWidth, lengthTarget);
199 gebp_kernel(res+startTarget, resStride, blockA, blockB, lengthTarget, actualPanelWidth, cols, alpha,
200 actualPanelWidth, actual_kc, 0, blockBOffset, blockW);
206 Index start = IsLower ? k2 : 0;
207 Index end = IsLower ? rows : (std::min)(actual_k2,rows);
208 for(Index i2=start; i2<end; i2+=mc)
210 const Index actual_mc = (std::min)(i2+mc,end)-i2;
212 (blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
214 gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1, -1, 0, 0, blockW);
221 template <
typename Scalar,
typename Index,
int Mode,
222 int LhsStorageOrder,
bool ConjugateLhs,
223 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
225 LhsStorageOrder,ConjugateLhs,
226 RhsStorageOrder,ConjugateRhs,ColMajor,Version>
236 Index _rows, Index _cols, Index _depth,
237 const Scalar* _lhs, Index lhsStride,
238 const Scalar* _rhs, Index rhsStride,
239 Scalar* res, Index resStride,
243 template <
typename Scalar,
typename Index,
int Mode,
244 int LhsStorageOrder,
bool ConjugateLhs,
245 int RhsStorageOrder,
bool ConjugateRhs,
int Version>
247 LhsStorageOrder,ConjugateLhs,
248 RhsStorageOrder,ConjugateRhs,
ColMajor,Version>::run(
249 Index _rows, Index _cols, Index _depth,
250 const Scalar* _lhs, Index lhsStride,
251 const Scalar* _rhs, Index rhsStride,
252 Scalar* res, Index resStride,
256 Index diagSize = (std::min)(_cols,_depth);
258 Index depth = IsLower ? _depth : diagSize;
259 Index cols = IsLower ? diagSize : _cols;
264 Index kc = blocking.
kc();
265 Index mc = (std::min)(rows,blocking.
mc());
267 std::size_t sizeA = kc*mc;
268 std::size_t sizeB = kc*cols;
269 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
278 triangularBuffer.diagonal().setZero();
280 triangularBuffer.diagonal().setOnes();
287 for(Index k2=IsLower ? 0 : depth;
288 IsLower ? k2<depth : k2>0;
289 IsLower ? k2+=kc : k2-=kc)
291 Index actual_kc = (std::min)(IsLower ? depth-k2 : k2, kc);
292 Index actual_k2 = IsLower ? k2 : k2-actual_kc;
295 if(IsLower && (k2<cols) && (actual_k2+actual_kc>cols))
298 k2 = actual_k2 + actual_kc - kc;
302 Index rs = IsLower ? (std::min)(cols,actual_k2) : cols - k2;
304 Index ts = (IsLower && actual_k2>=cols) ? 0 : actual_kc;
306 Scalar* geb = blockB+ts*ts;
308 pack_rhs(geb, &
rhs(actual_k2,IsLower ? 0 : k2), rhsStride, actual_kc, rs);
313 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
315 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
316 Index actual_j2 = actual_k2 + j2;
317 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
318 Index panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
320 pack_rhs_panel(blockB+j2*actual_kc,
321 &
rhs(actual_k2+panelOffset, actual_j2), rhsStride,
322 panelLength, actualPanelWidth,
323 actual_kc, panelOffset);
326 for (Index j=0;j<actualPanelWidth;++j)
329 triangularBuffer.coeffRef(j,j) =
rhs(actual_j2+j,actual_j2+j);
330 for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k)
331 triangularBuffer.coeffRef(k,j) =
rhs(actual_j2+k,actual_j2+j);
334 pack_rhs_panel(blockB+j2*actual_kc,
335 triangularBuffer.data(), triangularBuffer.outerStride(),
336 actualPanelWidth, actualPanelWidth,
341 for (Index i2=0; i2<rows; i2+=mc)
343 const Index actual_mc = (std::min)(mc,rows-i2);
344 pack_lhs(blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
349 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
351 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
352 Index panelLength = IsLower ? actual_kc-j2 : j2+actualPanelWidth;
353 Index blockOffset = IsLower ? j2 : 0;
355 gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride,
356 blockA, blockB+j2*actual_kc,
357 actual_mc, panelLength, actualPanelWidth,
359 actual_kc, actual_kc,
360 blockOffset, blockOffset,
364 gebp_kernel(res+i2+(IsLower ? 0 : k2)*resStride, resStride,
365 blockA, geb, actual_mc, actual_kc, rs,
367 -1, -1, 0, 0, blockW);
376 template<
int Mode,
bool LhsIsTriangular,
typename Lhs,
typename Rhs>
378 :
traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>, Lhs, Rhs> >
383 template<
int Mode,
bool LhsIsTriangular,
typename Lhs,
typename Rhs>
385 :
public ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>, Lhs, Rhs >
391 template<
typename Dest>
void scaleAndAddTo(Dest& dst,
const Scalar& alpha)
const 396 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
397 * RhsBlasTraits::extractScalarFactor(m_rhs);
400 Lhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxColsAtCompileTime,4> BlockingType;
403 Index stripedRows = ((!LhsIsTriangular) || (IsLower)) ? lhs.rows() : (std::min)(lhs.rows(),lhs.cols());
404 Index stripedCols = ((LhsIsTriangular) || (!IsLower)) ? rhs.cols() : (std::min)(rhs.cols(),rhs.rows());
405 Index stripedDepth = LhsIsTriangular ? ((!IsLower) ? lhs.cols() : (std::min)(lhs.cols(),lhs.rows()))
406 : ((IsLower) ? rhs.rows() : (std::min)(rhs.rows(),rhs.cols()));
408 BlockingType blocking(stripedRows, stripedCols, stripedDepth);
411 Mode, LhsIsTriangular,
416 stripedRows, stripedCols, stripedDepth,
417 &lhs.coeffRef(0,0), lhs.outerStride(),
418 &rhs.coeffRef(0,0), rhs.outerStride(),
419 &dst.coeffRef(0,0), dst.outerStride(),
420 actualAlpha, blocking
427 #endif // EIGEN_TRIANGULAR_MATRIX_MATRIX_H
#define EIGEN_PRODUCT_PUBLIC_INTERFACE(Derived)
#define EIGEN_STRONG_INLINE
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
gebp_traits< Scalar, Scalar > Traits
void scaleAndAddTo(Dest &dst, const Scalar &alpha) const
iterative scaling algorithm to equilibrate rows and column norms in matrices
gebp_traits< Scalar, Scalar > Traits
const unsigned int RowMajorBit
static EIGEN_STRONG_INLINE void run(Index rows, Index cols, Index depth, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsStride, Scalar *res, Index resStride, const Scalar &alpha, level3_blocking< Scalar, Scalar > &blocking)
Derived & setZero(Index size)
void rhs(const real_t *x, real_t *f)
#define EIGEN_DONT_INLINE
The matrix class, also used for vectors and row-vectors.
#define EIGEN_PLAIN_ENUM_MAX(a, b)