31 #ifndef SPARSELU_COLUMN_BMOD_H 32 #define SPARSELU_COLUMN_BMOD_H 52 template <
typename Scalar,
typename StorageIndex>
53 Index SparseLUImpl<Scalar,StorageIndex>::column_bmod(
const Index jcol,
const Index nseg, BlockScalarVector dense, ScalarVector& tempv,
54 BlockIndexVector segrep, BlockIndexVector repfnz,
Index fpanelc, GlobalLU_t& glu)
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;
93 kfnz = (std::max)(kfnz, fpanelc);
95 segsize = krep - kfnz + 1;
96 nsupc = krep - fst_col + 1;
97 nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc);
98 nrow = nsupr - d_fsupc - nsupc;
99 Index lda = glu.xlusup(fst_col+1) - glu.xlusup(fst_col);
104 no_zeros = kfnz - fst_col;
106 LU_kernel_bmod<1>::run(segsize, dense, tempv, glu.lusup, luptr, lda, nrow, glu.lsub, lptr, no_zeros);
108 LU_kernel_bmod<Dynamic>::run(segsize, dense, tempv, glu.lusup, luptr, lda, nrow, glu.lsub, lptr, no_zeros);
113 nextlu = glu.xlusup(jcol);
114 fsupc = glu.xsup(jsupno);
118 new_next = nextlu + glu.xlsub(fsupc + 1) - glu.xlsub(fsupc);
119 Index offset = internal::first_multiple<Index>(new_next, internal::packet_traits<Scalar>::size) - new_next;
122 while (new_next > glu.nzlumax )
124 mem = memXpand<ScalarVector>(glu.lusup, glu.nzlumax, nextlu, LUSUP, glu.num_expansions);
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);
138 glu.lusup.segment(nextlu,offset).setZero();
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;
165 Index lda = glu.xlusup(jcol+1) - glu.xlusup(jcol);
166 MappedMatrixBlock A( &(glu.lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
167 VectorBlock<ScalarVector> u(glu.lusup, ufirst, nsupc);
168 u = A.template triangularView<UnitLower>().solve(u);
170 new (&A) MappedMatrixBlock ( &(glu.lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
171 VectorBlock<ScalarVector> l(glu.lusup, ufirst+nsupc, nrow);
172 l.noalias() -= A * u;
181 #endif // SPARSELU_COLUMN_BMOD_H Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Eigen_Colamd.h:50