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)
41 triangular_solve_matrix<
46 ::run(size, cols, tri, triStride, _other, otherStride);
52 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
53 struct triangular_solve_matrix<Scalar,Index,
OnTheLeft,Mode,Conjugate,TriStorageOrder,
ColMajor>
56 Index size, Index otherSize,
57 const Scalar* _tri, Index triStride,
58 Scalar* _other, Index otherStride)
60 Index cols = otherSize;
61 const_blas_data_mapper<Scalar, Index, TriStorageOrder> tri(_tri,triStride);
62 blas_data_mapper<Scalar, Index, ColMajor> other(_other,otherStride);
64 typedef gebp_traits<Scalar,Scalar> Traits;
73 computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
75 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
76 std::size_t sizeB = sizeW + kc*cols;
79 Scalar* blockB = allocatedBlockB + sizeW;
80 Scalar* blockW = allocatedBlockB;
82 conj_if<Conjugate>
conj;
83 gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, Conjugate, false> gebp_kernel;
84 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, TriStorageOrder> pack_lhs;
85 gemm_pack_rhs<Scalar, Index, Traits::nr, ColMajor, false, true> pack_rhs;
89 std::ptrdiff_t l1, l2;
91 Index subcols = cols>0 ? l2/(4 *
sizeof(Scalar) * otherStride) : 0;
92 subcols = std::max<Index>((subcols/Traits::nr)*Traits::nr, Traits::nr);
94 for(Index k2=IsLower ? 0 : size;
95 IsLower ? k2<size : k2>0;
96 IsLower ? k2+=kc : k2-=kc)
98 const Index actual_kc = (std::min)(IsLower ? size-k2 : k2, kc);
113 for(Index j2=0; j2<cols; j2+=subcols)
115 Index actual_cols = (std::min)(cols-j2,subcols);
117 for (Index k1=0; k1<actual_kc; k1+=SmallPanelWidth)
119 Index actualPanelWidth = std::min<Index>(actual_kc-k1, SmallPanelWidth);
121 for (Index k=0; k<actualPanelWidth; ++k)
124 Index i = IsLower ? k2+k1+k : k2-k1-k-1;
125 Index s = IsLower ? k2+k1 : i+1;
126 Index rs = actualPanelWidth - k - 1;
128 Scalar a = (Mode &
UnitDiag) ? Scalar(1) : Scalar(1)/
conj(tri(i,i));
129 for (Index j=j2; j<j2+actual_cols; ++j)
134 const Scalar* l = &tri(i,s);
135 Scalar* r = &other(s,j);
136 for (Index i3=0; i3<k; ++i3)
137 b +=
conj(l[i3]) * r[i3];
139 other(i,j) = (other(i,j) - b)*a;
143 Index s = IsLower ? i+1 : i-rs;
144 Scalar b = (other(i,j) *= a);
145 Scalar* r = &other(s,j);
146 const Scalar* l = &tri(s,i);
147 for (Index i3=0;i3<rs;++i3)
148 r[i3] -= b *
conj(l[i3]);
153 Index lengthTarget = actual_kc-k1-actualPanelWidth;
154 Index startBlock = IsLower ? k2+k1 : k2-k1-actualPanelWidth;
155 Index blockBOffset = IsLower ? k1 : lengthTarget;
158 pack_rhs(blockB+actual_kc*j2, &other(startBlock,j2), otherStride, actualPanelWidth, actual_cols, actual_kc, blockBOffset);
163 Index startTarget = IsLower ? k2+k1+actualPanelWidth : k2-actual_kc;
165 pack_lhs(blockA, &tri(startTarget,startBlock), triStride, actualPanelWidth, lengthTarget);
167 gebp_kernel(&other(startTarget,j2), otherStride, blockA, blockB+actual_kc*j2, lengthTarget, actualPanelWidth, actual_cols, Scalar(-1),
168 actualPanelWidth, actual_kc, 0, blockBOffset, blockW);
175 Index start = IsLower ? k2+kc : 0;
176 Index end = IsLower ? size : k2-kc;
177 for(Index i2=start; i2<end; i2+=mc)
179 const Index actual_mc = (std::min)(mc,end-i2);
182 pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
184 gebp_kernel(_other+i2, otherStride, blockA, blockB, actual_mc, actual_kc, cols, Scalar(-1));
194 template <
typename Scalar,
typename Index,
int Mode,
bool Conjugate,
int TriStorageOrder>
195 struct triangular_solve_matrix<Scalar,Index,
OnTheRight,Mode,Conjugate,TriStorageOrder,
ColMajor>
198 Index size, Index otherSize,
199 const Scalar* _tri, Index triStride,
200 Scalar* _other, Index otherStride)
202 Index rows = otherSize;
203 const_blas_data_mapper<Scalar, Index, TriStorageOrder> rhs(_tri,triStride);
204 blas_data_mapper<Scalar, Index, ColMajor> lhs(_other,otherStride);
206 typedef gebp_traits<Scalar,Scalar> Traits;
208 RhsStorageOrder = TriStorageOrder,
219 computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
221 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
222 std::size_t sizeB = sizeW + kc*size;
225 Scalar* blockB = allocatedBlockB + sizeW;
227 conj_if<Conjugate>
conj;
228 gebp_kernel<Scalar,Scalar, Index, Traits::mr, Traits::nr, false, Conjugate> gebp_kernel;
229 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs;
230 gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder,false,true> pack_rhs_panel;
231 gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, ColMajor, false, true> pack_lhs_panel;
233 for(Index k2=IsLower ? size : 0;
234 IsLower ? k2>0 : k2<size;
235 IsLower ? k2-=kc : k2+=kc)
237 const Index actual_kc = (std::min)(IsLower ? k2 : size-k2, kc);
238 Index actual_k2 = IsLower ? k2-actual_kc : k2 ;
240 Index startPanel = IsLower ? 0 : k2+actual_kc;
241 Index rs = IsLower ? actual_k2 : size - actual_k2 - actual_kc;
242 Scalar* geb = blockB+actual_kc*actual_kc;
244 if (rs>0) pack_rhs(geb, &rhs(actual_k2,startPanel), triStride, actual_kc, rs);
249 for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
251 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
252 Index actual_j2 = actual_k2 + j2;
253 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
254 Index panelLength = IsLower ? actual_kc-j2-actualPanelWidth : j2;
257 pack_rhs_panel(blockB+j2*actual_kc,
258 &rhs(actual_k2+panelOffset, actual_j2), triStride,
259 panelLength, actualPanelWidth,
260 actual_kc, panelOffset);
264 for(Index i2=0; i2<rows; i2+=mc)
266 const Index actual_mc = (std::min)(mc,rows-i2);
271 for (Index j2 = IsLower
272 ? (actual_kc - ((actual_kc%SmallPanelWidth) ? Index(actual_kc%SmallPanelWidth)
273 : Index(SmallPanelWidth)))
275 IsLower ? j2>=0 : j2<actual_kc;
276 IsLower ? j2-=SmallPanelWidth : j2+=SmallPanelWidth)
278 Index actualPanelWidth = std::min<Index>(actual_kc-j2, SmallPanelWidth);
279 Index absolute_j2 = actual_k2 + j2;
280 Index panelOffset = IsLower ? j2+actualPanelWidth : 0;
281 Index panelLength = IsLower ? actual_kc - j2 - actualPanelWidth : j2;
286 gebp_kernel(&lhs(i2,absolute_j2), otherStride,
287 blockA, blockB+j2*actual_kc,
288 actual_mc, panelLength, actualPanelWidth,
290 actual_kc, actual_kc,
291 panelOffset, panelOffset,
296 for (Index k=0; k<actualPanelWidth; ++k)
298 Index j = IsLower ? absolute_j2+actualPanelWidth-k-1 : absolute_j2+k;
300 Scalar* r = &lhs(i2,j);
301 for (Index k3=0; k3<k; ++k3)
303 Scalar b =
conj(rhs(IsLower ? j+1+k3 : absolute_j2+k3,j));
304 Scalar* a = &lhs(i2,IsLower ? j+1+k3 : absolute_j2+k3);
305 for (Index i=0; i<actual_mc; ++i)
308 Scalar b = (Mode &
UnitDiag) ? Scalar(1) : Scalar(1)/
conj(rhs(j,j));
309 for (Index i=0; i<actual_mc; ++i)
314 pack_lhs_panel(blockA, _other+absolute_j2*otherStride+i2, otherStride,
315 actualPanelWidth, actual_mc,
321 gebp_kernel(_other+i2+startPanel*otherStride, otherStride, blockA, geb,
322 actual_mc, actual_kc, rs, Scalar(-1),
323 -1, -1, 0, 0, allocatedBlockB);
333 #endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_H