25 #ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_H
26 #define EIGEN_TRIANGULAR_SOLVER_MATRIX_H
33 template <
typename Scalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int TriStorageOrder>
34 struct triangular_solve_matrix<Scalar,Index,Side,Mode,Conjugate,TriStorageOrder,
RowMajor>
37 Index size, Index cols,
38 const Scalar* tri, Index triStride,
39 Scalar* _other, Index otherStride,
40 level3_blocking<Scalar,Scalar>& blocking)
42 triangular_solve_matrix<
47 ::run(size, cols, tri, triStride, _other, otherStride, blocking);
53 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
54 struct triangular_solve_matrix<Scalar,Index,
OnTheLeft,Mode,Conjugate,TriStorageOrder,
ColMajor>
57 Index size, Index otherSize,
58 const Scalar* _tri, Index triStride,
59 Scalar* _other, Index otherStride,
60 level3_blocking<Scalar,Scalar>& blocking)
62 Index cols = otherSize;
63 const_blas_data_mapper<Scalar, Index, TriStorageOrder> tri(_tri,triStride);
64 blas_data_mapper<Scalar, Index, ColMajor> other(_other,otherStride);
66 typedef gebp_traits<Scalar,Scalar> Traits;
72 Index kc = blocking.kc();
73 Index mc = (std::min)(size,blocking.mc());
75 std::size_t sizeA = kc*mc;
76 std::size_t sizeB = kc*cols;
77 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
83 conj_if<Conjugate>
conj;
84 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, Conjugate, false> gebp_kernel;
85 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, TriStorageOrder> pack_lhs;
86 gemm_pack_rhs<Scalar, Index, Traits::nr, ColMajor, false, true> pack_rhs;
90 std::ptrdiff_t l1, l2;
92 Index subcols = cols>0 ? l2/(4 *
sizeof(Scalar) * otherStride) : 0;
93 subcols = std::max<Index>((subcols/Traits::nr)*Traits::nr, Traits::nr);
95 for(Index k2=IsLower ? 0 : size;
96 IsLower ? k2<size : k2>0;
97 IsLower ? k2+=kc : k2-=kc)
99 const Index actual_kc = (std::min)(IsLower ? size-k2 : k2, kc);
114 for(Index j2=0; j2<cols; j2+=subcols)
116 Index actual_cols = (std::min)(cols-j2,subcols);
118 for (Index k1=0; k1<actual_kc; k1+=SmallPanelWidth)
120 Index actualPanelWidth = std::min<Index>(actual_kc-k1, SmallPanelWidth);
122 for (Index k=0; k<actualPanelWidth; ++k)
125 Index i = IsLower ? k2+k1+k : k2-k1-k-1;
126 Index s = IsLower ? k2+k1 : i+1;
127 Index rs = actualPanelWidth - k - 1;
129 Scalar a = (Mode &
UnitDiag) ? Scalar(1) : Scalar(1)/
conj(tri(i,i));
130 for (Index j=j2; j<j2+actual_cols; ++j)
135 const Scalar* l = &tri(i,s);
136 Scalar* r = &other(s,j);
137 for (Index i3=0; i3<k; ++i3)
138 b +=
conj(l[i3]) * r[i3];
140 other(i,j) = (other(i,j) - b)*a;
144 Index s = IsLower ? i+1 : i-rs;
145 Scalar b = (other(i,j) *= a);
146 Scalar* r = &other(s,j);
147 const Scalar* l = &tri(s,i);
148 for (Index i3=0;i3<rs;++i3)
149 r[i3] -= b *
conj(l[i3]);
154 Index lengthTarget = actual_kc-k1-actualPanelWidth;
155 Index startBlock = IsLower ? k2+k1 : k2-k1-actualPanelWidth;
156 Index blockBOffset = IsLower ? k1 : lengthTarget;
159 pack_rhs(blockB+actual_kc*j2, &other(startBlock,j2), otherStride, actualPanelWidth, actual_cols, actual_kc, blockBOffset);
164 Index startTarget = IsLower ? k2+k1+actualPanelWidth : k2-actual_kc;
166 pack_lhs(blockA, &tri(startTarget,startBlock), triStride, actualPanelWidth, lengthTarget);
168 gebp_kernel(&other(startTarget,j2), otherStride, blockA, blockB+actual_kc*j2, lengthTarget, actualPanelWidth, actual_cols, Scalar(-1),
169 actualPanelWidth, actual_kc, 0, blockBOffset, blockW);
176 Index start = IsLower ? k2+kc : 0;
177 Index end = IsLower ? size : k2-kc;
178 for(Index i2=start; i2<end; i2+=mc)
180 const Index actual_mc = (std::min)(mc,end-i2);
183 pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
185 gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1), -1, -1, 0, 0, blockW);
195 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
196 struct triangular_solve_matrix<Scalar,Index,
OnTheRight,Mode,Conjugate,TriStorageOrder,
ColMajor>
199 Index size, Index otherSize,
200 const Scalar* _tri, Index triStride,
201 Scalar* _other, Index otherStride,
202 level3_blocking<Scalar,Scalar>& blocking)
204 Index rows = otherSize;
205 const_blas_data_mapper<Scalar, Index, TriStorageOrder> rhs(_tri,triStride);
206 blas_data_mapper<Scalar, Index, ColMajor> lhs(_other,otherStride);
208 typedef gebp_traits<Scalar,Scalar> Traits;
210 RhsStorageOrder = TriStorageOrder,
215 Index kc = blocking.kc();
216 Index mc = (std::min)(rows,blocking.mc());
218 std::size_t sizeA = kc*mc;
219 std::size_t sizeB = kc*size;
220 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
226 conj_if<Conjugate>
conj;
227 gebp_kernel<Scalar,Scalar, Index, Traits::mr, Traits::nr, false, Conjugate> gebp_kernel;
228 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
229 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder,false,true> pack_rhs_panel;
230 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, ColMajor, false, true> pack_lhs_panel;
232 for(Index k2=IsLower ? size : 0;
233 IsLower ? k2>0 : k2<size;
234 IsLower ? k2-=kc : k2+=kc)
236 const Index actual_kc = (std::min)(IsLower ? k2 : size-k2, kc);
237 Index actual_k2 = IsLower ? k2-actual_kc : k2 ;
239 Index startPanel = IsLower ? 0 : k2+actual_kc;
240 Index rs = IsLower ? actual_k2 : size - actual_k2 - actual_kc;
241 Scalar* geb = blockB+actual_kc*actual_kc;
243 if (rs>0) pack_rhs(geb, &rhs(actual_k2,startPanel), triStride, actual_kc, rs);
248 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
250 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
251 Index actual_j2 = actual_k2 + j2;
252 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
253 Index panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
256 pack_rhs_panel(blockB+j2*actual_kc,
257 &rhs(actual_k2+panelOffset, actual_j2), triStride,
258 panelLength, actualPanelWidth,
259 actual_kc, panelOffset);
263 for(Index i2=0; i2<rows; i2+=mc)
265 const Index actual_mc = (std::min)(mc,rows-i2);
270 for (Index j2 = IsLower
271 ? (actual_kc - ((actual_kc%SmallPanelWidth) ? Index(actual_kc%SmallPanelWidth)
272 : Index(SmallPanelWidth)))
274 IsLower ? j2>=0 : j2<actual_kc;
275 IsLower ? j2-=SmallPanelWidth : j2+=SmallPanelWidth)
277 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
278 Index absolute_j2 = actual_k2 + j2;
279 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
280 Index panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
285 gebp_kernel(&lhs(i2,absolute_j2), otherStride,
286 blockA, blockB+j2*actual_kc,
287 actual_mc, panelLength, actualPanelWidth,
289 actual_kc, actual_kc,
290 panelOffset, panelOffset,
295 for (Index k=0; k<actualPanelWidth; ++k)
297 Index j = IsLower ? absolute_j2+actualPanelWidth-k-1 : absolute_j2+k;
299 Scalar* r = &lhs(i2,j);
300 for (Index k3=0; k3<k; ++k3)
302 Scalar b =
conj(rhs(IsLower ? j+1+k3 : absolute_j2+k3,j));
303 Scalar* a = &lhs(i2,IsLower ? j+1+k3 : absolute_j2+k3);
304 for (Index i=0; i<actual_mc; ++i)
307 Scalar b = (Mode &
UnitDiag) ? Scalar(1) : Scalar(1)/
conj(rhs(j,j));
308 for (Index i=0; i<actual_mc; ++i)
313 pack_lhs_panel(blockA, _other+absolute_j2*otherStride+i2, otherStride,
314 actualPanelWidth, actual_mc,
320 gebp_kernel(_other+i2+startPanel*otherStride, otherStride, blockA, geb,
321 actual_mc, actual_kc, rs, Scalar(-1),
322 -1, -1, 0, 0, blockW);
332 #endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_H