10 #ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_H 11 #define EIGEN_TRIANGULAR_SOLVER_MATRIX_H 18 template <
typename Scalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int TriStorageOrder>
22 Index size, Index cols,
23 const Scalar* tri, Index triStride,
24 Scalar* _other, Index otherStride,
32 ::run(size, cols, tri, triStride, _other, otherStride, blocking);
38 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
42 Index size, Index otherSize,
43 const Scalar* _tri, Index triStride,
44 Scalar* _other, Index otherStride,
47 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
49 Index size, Index otherSize,
50 const Scalar* _tri, Index triStride,
51 Scalar* _other, Index otherStride,
54 Index cols = otherSize;
64 Index kc = blocking.
kc();
65 Index mc = (std::min)(size,blocking.
mc());
67 std::size_t sizeA = kc*mc;
68 std::size_t sizeB = kc*cols;
69 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
82 std::ptrdiff_t l1, l2;
84 Index subcols = cols>0 ? l2/(4 *
sizeof(Scalar) * otherStride) : 0;
85 subcols = std::max<Index>((subcols/Traits::nr)*Traits::nr, Traits::nr);
87 for(Index k2=IsLower ? 0 : size;
88 IsLower ? k2<size : k2>0;
89 IsLower ? k2+=kc : k2-=kc)
91 const Index actual_kc = (std::min)(IsLower ? size-k2 : k2, kc);
106 for(Index j2=0; j2<cols; j2+=subcols)
108 Index actual_cols = (std::min)(cols-j2,subcols);
110 for (Index k1=0; k1<actual_kc; k1+=SmallPanelWidth)
112 Index actualPanelWidth = std::min<Index>(actual_kc-k1, SmallPanelWidth);
114 for (Index k=0; k<actualPanelWidth; ++k)
117 Index i = IsLower ? k2+k1+k : k2-k1-k-1;
118 Index s = IsLower ? k2+k1 : i+1;
119 Index rs = actualPanelWidth - k - 1;
121 Scalar a = (Mode &
UnitDiag) ? Scalar(1) : Scalar(1)/
conj(tri(i,i));
122 for (Index j=j2; j<j2+actual_cols; ++j)
127 const Scalar* l = &tri(i,s);
128 Scalar* r = &other(s,j);
129 for (Index i3=0; i3<k; ++i3)
130 b +=
conj(l[i3]) * r[i3];
132 other(i,j) = (other(i,j) - b)*a;
136 Index s = IsLower ? i+1 : i-rs;
137 Scalar b = (other(i,j) *= a);
138 Scalar* r = &other(s,j);
139 const Scalar* l = &tri(s,i);
140 for (Index i3=0;i3<rs;++i3)
141 r[i3] -= b *
conj(l[i3]);
146 Index lengthTarget = actual_kc-k1-actualPanelWidth;
147 Index startBlock = IsLower ? k2+k1 : k2-k1-actualPanelWidth;
148 Index blockBOffset = IsLower ? k1 : lengthTarget;
151 pack_rhs(blockB+actual_kc*j2, &other(startBlock,j2), otherStride, actualPanelWidth, actual_cols, actual_kc, blockBOffset);
156 Index startTarget = IsLower ? k2+k1+actualPanelWidth : k2-actual_kc;
158 pack_lhs(blockA, &tri(startTarget,startBlock), triStride, actualPanelWidth, lengthTarget);
160 gebp_kernel(&other(startTarget,j2), otherStride, blockA, blockB+actual_kc*j2, lengthTarget, actualPanelWidth, actual_cols, Scalar(-1),
161 actualPanelWidth, actual_kc, 0, blockBOffset, blockW);
168 Index start = IsLower ? k2+kc : 0;
169 Index end = IsLower ? size : k2-kc;
170 for(Index i2=start; i2<end; i2+=mc)
172 const Index actual_mc = (std::min)(mc,end-i2);
175 pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
177 gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1), -1, -1, 0, 0, blockW);
186 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
190 Index size, Index otherSize,
191 const Scalar* _tri, Index triStride,
192 Scalar* _other, Index otherStride,
195 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
197 Index size, Index otherSize,
198 const Scalar* _tri, Index triStride,
199 Scalar* _other, Index otherStride,
202 Index rows = otherSize;
208 RhsStorageOrder = TriStorageOrder,
213 Index kc = blocking.
kc();
214 Index mc = (std::min)(rows,blocking.
mc());
216 std::size_t sizeA = kc*mc;
217 std::size_t sizeB = kc*size;
218 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
230 for(Index k2=IsLower ? size : 0;
231 IsLower ? k2>0 : k2<size;
232 IsLower ? k2-=kc : k2+=kc)
234 const Index actual_kc = (std::min)(IsLower ? k2 : size-k2, kc);
235 Index actual_k2 = IsLower ? k2-actual_kc : k2 ;
237 Index startPanel = IsLower ? 0 : k2+actual_kc;
238 Index rs = IsLower ? actual_k2 : size - actual_k2 - actual_kc;
239 Scalar* geb = blockB+actual_kc*actual_kc;
241 if (rs>0) pack_rhs(geb, &
rhs(actual_k2,startPanel), triStride, actual_kc, rs);
246 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
248 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
249 Index actual_j2 = actual_k2 + j2;
250 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
251 Index panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
254 pack_rhs_panel(blockB+j2*actual_kc,
255 &
rhs(actual_k2+panelOffset, actual_j2), triStride,
256 panelLength, actualPanelWidth,
257 actual_kc, panelOffset);
261 for(Index i2=0; i2<rows; i2+=mc)
263 const Index actual_mc = (std::min)(mc,rows-i2);
268 for (Index j2 = IsLower
269 ? (actual_kc - ((actual_kc%SmallPanelWidth) ? Index(actual_kc%SmallPanelWidth)
270 : Index(SmallPanelWidth)))
272 IsLower ? j2>=0 : j2<actual_kc;
273 IsLower ? j2-=SmallPanelWidth : j2+=SmallPanelWidth)
275 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
276 Index absolute_j2 = actual_k2 + j2;
277 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
278 Index panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
283 gebp_kernel(&lhs(i2,absolute_j2), otherStride,
284 blockA, blockB+j2*actual_kc,
285 actual_mc, panelLength, actualPanelWidth,
287 actual_kc, actual_kc,
288 panelOffset, panelOffset,
293 for (Index k=0; k<actualPanelWidth; ++k)
295 Index j = IsLower ? absolute_j2+actualPanelWidth-k-1 : absolute_j2+k;
297 Scalar* r = &lhs(i2,j);
298 for (Index k3=0; k3<k; ++k3)
300 Scalar b =
conj(
rhs(IsLower ? j+1+k3 : absolute_j2+k3,j));
301 Scalar* a = &lhs(i2,IsLower ? j+1+k3 : absolute_j2+k3);
302 for (Index i=0; i<actual_mc; ++i)
306 for (Index i=0; i<actual_mc; ++i)
311 pack_lhs_panel(blockA, _other+absolute_j2*otherStride+i2, otherStride,
312 actualPanelWidth, actual_mc,
318 gebp_kernel(_other+i2+startPanel*otherStride, otherStride, blockA, geb,
319 actual_mc, actual_kc, rs, Scalar(-1),
320 -1, -1, 0, 0, blockW);
329 #endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_H
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
iterative scaling algorithm to equilibrate rows and column norms in matrices
void manage_caching_sizes(Action action, std::ptrdiff_t *l1=0, std::ptrdiff_t *l2=0)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
static void run(Index size, Index cols, const Scalar *tri, Index triStride, Scalar *_other, Index otherStride, level3_blocking< Scalar, Scalar > &blocking)
void rhs(const real_t *x, real_t *f)
#define EIGEN_DONT_INLINE
#define EIGEN_PLAIN_ENUM_MAX(a, b)