31 #ifndef SPARSELU_COLUMN_BMOD_H
32 #define SPARSELU_COLUMN_BMOD_H
52 template <
typename Scalar,
typename StorageIndex>
56 Index jsupno, k, ksub, krep, ksupno;
57 Index lptr, nrow, isub, irow, nextlu, new_next, ufirst;
58 Index fsupc, nsupc, nsupr, luptr, kfnz, no_zeros;
68 jsupno = glu.
supno(jcol);
75 for (ksub = 0; ksub < nseg; ksub++)
77 krep = segrep(k); k--;
78 ksupno = glu.
supno(krep);
79 if (jsupno != ksupno )
82 fsupc = glu.
xsup(ksupno);
83 fst_col = (
std::max)(fsupc, fpanelc);
87 d_fsupc = fst_col - fsupc;
89 luptr = glu.
xlusup(fst_col) + d_fsupc;
90 lptr = glu.
xlsub(fsupc) + d_fsupc;
95 segsize = krep - kfnz + 1;
96 nsupc = krep - fst_col + 1;
98 nrow = nsupr - d_fsupc - nsupc;
104 no_zeros = kfnz - fst_col;
113 nextlu = glu.
xlusup(jcol);
114 fsupc = glu.
xsup(jsupno);
118 new_next = nextlu + glu.
xlsub(fsupc + 1) - glu.
xlsub(fsupc);
122 while (new_next > glu.
nzlumax )
128 for (isub = glu.
xlsub(fsupc); isub < glu.
xlsub(fsupc+1); isub++)
130 irow = glu.
lsub(isub);
131 glu.
lusup(nextlu) = dense(irow);
132 dense(irow) =
Scalar(0.0);
141 glu.
xlusup(jcol + 1) = StorageIndex(nextlu);
149 fst_col = (
std::max)(fsupc, fpanelc);
155 d_fsupc = fst_col - fsupc;
157 lptr = glu.
xlsub(fsupc) + d_fsupc;
158 luptr = glu.
xlusup(fst_col) + d_fsupc;
159 nsupr = glu.
xlsub(fsupc+1) - glu.
xlsub(fsupc);
160 nsupc = jcol - fst_col;
161 nrow = nsupr - d_fsupc - nsupc;
164 ufirst = glu.
xlusup(jcol) + d_fsupc;
168 u =
A.template triangularView<UnitLower>().solve(u);
172 l.noalias() -=
A * u;
181 #endif // SPARSELU_COLUMN_BMOD_H