30#ifndef SPARSELU_PANEL_DFS_H
31#define SPARSELU_PANEL_DFS_H
37template<
typename IndexVector>
38struct panel_dfs_traits
40 typedef typename IndexVector::Scalar StorageIndex;
41 panel_dfs_traits(
Index jcol, StorageIndex* marker)
42 : m_jcol(jcol), m_marker(marker)
44 bool update_segrep(
Index krep, StorageIndex jj)
46 if(m_marker[krep]<m_jcol)
53 void mem_expand(IndexVector& ,
Index ,
Index ) {}
54 enum { ExpandMem =
false };
56 StorageIndex* m_marker;
60template <
typename Scalar,
typename StorageIndex>
61template <
typename Traits>
62void SparseLUImpl<Scalar,StorageIndex>::dfs_kernel(
const StorageIndex jj, IndexVector& perm_r,
63 Index& nseg, IndexVector& panel_lsub, IndexVector& segrep,
64 Ref<IndexVector> repfnz_col, IndexVector& xprune, Ref<IndexVector> marker, IndexVector& parent,
65 IndexVector& xplore, GlobalLU_t& glu,
70 StorageIndex kmark = marker(krow);
74 StorageIndex kperm = perm_r(krow);
75 if (kperm == emptyIdxLU ) {
77 panel_lsub(nextl_col++) = StorageIndex(krow);
79 traits.mem_expand(panel_lsub, nextl_col, kmark);
86 StorageIndex krep = glu.xsup(glu.supno(kperm)+1) - 1;
88 StorageIndex myfnz = repfnz_col(krep);
90 if (myfnz != emptyIdxLU )
93 if (myfnz > kperm ) repfnz_col(krep) = kperm;
99 StorageIndex oldrep = emptyIdxLU;
100 parent(krep) = oldrep;
101 repfnz_col(krep) = kperm;
102 StorageIndex xdfs = glu.xlsub(krep);
103 Index maxdfs = xprune(krep);
109 while (xdfs < maxdfs)
111 StorageIndex kchild = glu.lsub(xdfs);
113 StorageIndex chmark = marker(kchild);
118 StorageIndex chperm = perm_r(kchild);
120 if (chperm == emptyIdxLU)
123 panel_lsub(nextl_col++) = kchild;
124 traits.mem_expand(panel_lsub, nextl_col, chmark);
131 StorageIndex chrep = glu.xsup(glu.supno(chperm)+1) - 1;
132 myfnz = repfnz_col(chrep);
134 if (myfnz != emptyIdxLU)
137 repfnz_col(chrep) = chperm;
144 parent(krep) = oldrep;
145 repfnz_col(krep) = chperm;
146 xdfs = glu.xlsub(krep);
147 maxdfs = xprune(krep);
160 if(traits.update_segrep(krep,jj))
169 if (kpar == emptyIdxLU)
173 maxdfs = xprune(krep);
175 }
while (kpar != emptyIdxLU);
218template <
typename Scalar,
typename StorageIndex>
219void SparseLUImpl<Scalar,StorageIndex>::panel_dfs(
const Index m,
const Index w,
const Index jcol, MatrixType& A, IndexVector& perm_r,
Index& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
224 VectorBlock<IndexVector> marker1(marker, m, m);
227 panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
230 for (StorageIndex jj = StorageIndex(jcol); jj < jcol + w; jj++)
232 nextl_col = (jj - jcol) * m;
234 VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m);
235 VectorBlock<ScalarVector> dense_col(dense,nextl_col, m);
239 for (
typename MatrixType::InnerIterator it(A, jj); it; ++it)
241 Index krow = it.row();
242 dense_col(krow) = it.value();
244 StorageIndex kmark = marker(krow);
248 dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent,
249 xplore, glu, nextl_col, krow, traits);
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