diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h index 36c97f947..a4562af9c 100644 --- a/Eigen/src/SparseLU/SparseLU_column_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h @@ -56,157 +56,68 @@ * > 0 number of bytes allocated when run out of space * */ +template +struct LU_column_dfs_traits +{ + typedef typename IndexVector::Scalar Index; + LU_column_dfs_traits(Index jcol, Index& jsuper, LU_GlobalLU_t& glu) + : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu) + {} + bool update_segrep(Index /*krep*/, Index /*jj*/) + { + return true; + } + void mem_expand(IndexVector& lsub, int& nextl, int chmark) + { + if (nextl >= m_glu.nzlmax) + LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); + if (chmark != (m_jcol-1)) m_jsuper_ref = IND_EMPTY; + } + enum { ExpandMem = true }; + + int m_jcol; + int& m_jsuper_ref; + LU_GlobalLU_t& m_glu; +}; + template int LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg, BlockIndexVector& lsub_col, IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu) { typedef typename IndexVector::Scalar Index; typedef typename ScalarVector::Scalar Scalar; - int jsuper, nsuper, nextl; - int krow; // Row index of the current element - int kperm; // permuted row index - int krep; // Supernode reprentative of the current row - int k, kmark; - int chperm, chmark, chrep, oldrep, kchild; - int myfnz; // First nonzero element in the current column - int xdfs, maxdfs, kpar; - int mem; // Initialize pointers IndexVector& xsup = glu.xsup; IndexVector& supno = glu.supno; IndexVector& lsub = glu.lsub; - IndexVector& xlsub = glu.xlsub; - Index& nzlmax = glu.nzlmax; + IndexVector& xlsub = glu.xlsub; - int jcolm1 = jcol - 1; - int jcolp1 = jcol + 1; - nsuper = supno(jcol); - jsuper = nsuper; - nextl = xlsub(jcol); + int jsuper = supno(jcol); + int nextl = xlsub(jcol); VectorBlock marker2(marker, 2*m, m); - int fsupc, jptr, jm1ptr, ito, ifrom, istop; + + + LU_column_dfs_traits traits(jcol, jsuper, glu); + // For each nonzero in A(*,jcol) do dfs - for (k = 0; lsub_col[k] != IND_EMPTY; k++) + for (int k = 0; lsub_col[k] != IND_EMPTY; k++) { - krow = lsub_col(k); + int krow = lsub_col(k); lsub_col(k) = IND_EMPTY; - kmark = marker2(krow); + int kmark = marker2(krow); // krow was visited before, go to the next nonz; - if (kmark == jcol) continue; - - // For each unmarker nbr krow of jcol - marker2(krow) = jcol; - kperm = perm_r(krow); - - if (kperm == IND_EMPTY ) - { - // krow is in L: place it in structure of L(*,jcol) - lsub(nextl++) = krow; // krow is indexed into A - if ( nextl >= nzlmax ) - { - mem = LUMemXpand(lsub, nzlmax, nextl, LSUB, glu.num_expansions); - if ( mem ) return mem; - } - if (kmark != jcolm1) jsuper = IND_EMPTY; // Row index subset testing - } - else - { - // krow is in U : if its supernode-rep krep - // has been explored, update repfnz(*) - krep = xsup(supno(kperm)+1) - 1; - myfnz = repfnz(krep); - - if (myfnz != IND_EMPTY ) - { - // visited before - if (myfnz > kperm) repfnz(krep) = kperm; - // continue; - } - else - { - // otherwise, perform dfs starting at krep - oldrep = IND_EMPTY; - parent(krep) = oldrep; - repfnz(krep) = kperm; - xdfs = xlsub(krep); - maxdfs = xprune(krep); - - do - { - // For each unmarked kchild of krep - while (xdfs < maxdfs) - { - kchild = lsub(xdfs); - xdfs++; - chmark = marker2(kchild); - - if (chmark != jcol) - { - // Not reached yet - marker2(kchild) = jcol; - chperm = perm_r(kchild); - - if (chperm == IND_EMPTY) - { - // if kchild is in L: place it in L(*,k) - lsub(nextl++) = kchild; - if (nextl >= nzlmax) - { - mem = LUMemXpand(lsub, nzlmax, nextl, LSUB, glu.num_expansions); - if (mem) return mem; - } - if (chmark != jcolm1) jsuper = IND_EMPTY; - } - else - { - // if kchild is in U : - // chrep = its supernode-rep. If its rep has been explored, - // update its repfnz - chrep = xsup(supno(chperm)+1) - 1; - myfnz = repfnz(chrep); - if (myfnz != IND_EMPTY) - { - // Visited before - if ( myfnz > chperm) repfnz(chrep) = chperm; - } - else - { - // continue dfs at super-rep of kchild - xplore(krep) = xdfs; - oldrep = krep; - krep = chrep; // Go deeped down G(L^t) - parent(krep) = oldrep; - repfnz(krep) = chperm; - xdfs = xlsub(krep); - maxdfs = xprune(krep); - } // else myfnz - } // else for chperm - - } // if chmark - - } // end while - - // krow has no more unexplored nbrs; - // place supernode-rep krep in postorder DFS. - // backtrack dfs to its parent - - segrep(nseg) = krep; - ++nseg; - kpar = parent(krep); // Pop from stack, mimic recursion - if (kpar == IND_EMPTY) break; // dfs done - krep = kpar; - xdfs = xplore(krep); - maxdfs = xprune(krep); - - } while ( kpar != IND_EMPTY); - - } // else myfnz - - } // else kperm + if (kmark == jcol) continue; + LU_dfs_kernel(jcol, perm_r, nseg, lsub, segrep, repfnz, xprune, marker2, parent, + xplore, glu, nextl, krow, traits); } // for each nonzero ... + int fsupc, jptr, jm1ptr, ito, ifrom, istop; + int nsuper = supno(jcol); + int jcolp1 = jcol + 1; + int jcolm1 = jcol - 1; + // check to see if j belongs in the same supernode as j-1 if ( jcol == 0 ) { // Do nothing for column 0 diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 79dd4da40..75fbd0b0e 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -29,6 +29,132 @@ */ #ifndef SPARSELU_PANEL_DFS_H #define SPARSELU_PANEL_DFS_H + +template +void LU_dfs_kernel(const int jj, IndexVector& perm_r, + int& nseg, IndexVector& panel_lsub, IndexVector& segrep, + VectorBlock& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent, + IndexVector& xplore, LU_GlobalLU_t& glu, + int& nextl_col, int krow, Traits& traits + ) +{ + IndexVector& xsup = glu.xsup; + IndexVector& supno = glu.supno; + IndexVector& lsub = glu.lsub; + IndexVector& xlsub = glu.xlsub; + + int kmark = marker(krow); + + // For each unmarked krow of jj + marker(krow) = jj; + int kperm = perm_r(krow); + if (kperm == IND_EMPTY ) { + // krow is in L : place it in structure of L(*, jj) + panel_lsub(nextl_col++) = krow; // krow is indexed into A + + traits.mem_expand(panel_lsub, nextl_col, kmark); + } + else + { + // krow is in U : if its supernode-representative krep + // has been explored, update repfnz(*) + // krep = supernode representative of the current row + int krep = xsup(supno(kperm)+1) - 1; + // First nonzero element in the current column: + int myfnz = repfnz_col(krep); + + if (myfnz != IND_EMPTY ) + { + // Representative visited before + if (myfnz > kperm ) repfnz_col(krep) = kperm; + + } + else + { + // Otherwise, perform dfs starting at krep + int oldrep = IND_EMPTY; + parent(krep) = oldrep; + repfnz_col(krep) = kperm; + int xdfs = xlsub(krep); + int maxdfs = xprune(krep); + + int kpar; + do + { + // For each unmarked kchild of krep + while (xdfs < maxdfs) + { + int kchild = lsub(xdfs); + xdfs++; + int chmark = marker(kchild); + + if (chmark != jj ) + { + marker(kchild) = jj; + int chperm = perm_r(kchild); + + if (chperm == IND_EMPTY) + { + // case kchild is in L: place it in L(*, j) + panel_lsub(nextl_col++) = kchild; + traits.mem_expand(panel_lsub, nextl_col, chmark); + } + else + { + // case kchild is in U : + // chrep = its supernode-rep. If its rep has been explored, + // update its repfnz(*) + int chrep = xsup(supno(chperm)+1) - 1; + myfnz = repfnz_col(chrep); + + if (myfnz != IND_EMPTY) + { // Visited before + if (myfnz > chperm) + repfnz_col(chrep) = chperm; + } + else + { // Cont. dfs at snode-rep of kchild + xplore(krep) = xdfs; + oldrep = krep; + krep = chrep; // Go deeper down G(L) + parent(krep) = oldrep; + repfnz_col(krep) = chperm; + xdfs = xlsub(krep); + maxdfs = xprune(krep); + + } // end if myfnz != -1 + } // end if chperm == -1 + + } // end if chmark !=jj + } // end while xdfs < maxdfs + + // krow has no more unexplored nbrs : + // Place snode-rep krep in postorder DFS, if this + // segment is seen for the first time. (Note that + // "repfnz(krep)" may change later.) + // Baktrack dfs to its parent + if(traits.update_segrep(krep,jj)) + //if (marker1(krep) < jcol ) + { + segrep(nseg) = krep; + ++nseg; + //marker1(krep) = jj; + } + + kpar = parent(krep); // Pop recursion, mimic recursion + if (kpar == IND_EMPTY) + break; // dfs done + krep = kpar; + xdfs = xplore(krep); + maxdfs = xprune(krep); + + } while (kpar != IND_EMPTY); // Do until empty stack + + } // end if (myfnz = -1) + + } // end if (kperm == -1) +} + /** * \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w) * @@ -62,29 +188,42 @@ * * */ + +template +struct LU_panel_dfs_traits +{ + typedef typename IndexVector::Scalar Index; + LU_panel_dfs_traits(Index jcol, Index* marker) + : m_jcol(jcol), m_marker(marker) + {} + bool update_segrep(Index krep, Index jj) + { + if(m_marker[krep] void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu) { - - int jj; // Index through each column in the panel int nextl_col; // Next available position in panel_lsub[*,jj] - int krow; // Row index of the current element - int kperm; // permuted row index - int krep; // Supernode representative of the current row - int kmark; - int chperm, chmark, chrep, oldrep, kchild; - int myfnz; // First nonzero element in the current column - int xdfs, maxdfs, kpar; // Initialize pointers VectorBlock marker1(marker, m, m); nseg = 0; - IndexVector& xsup = glu.xsup; - IndexVector& supno = glu.supno; - IndexVector& lsub = glu.lsub; - IndexVector& xlsub = glu.xlsub; + + LU_panel_dfs_traits traits(jcol, marker1.data()); + // For each column in the panel - for (jj = jcol; jj < jcol + w; jj++) + for (int jj = jcol; jj < jcol + w; jj++) { nextl_col = (jj - jcol) * m; @@ -95,114 +234,15 @@ void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, Index // For each nnz in A[*, jj] do depth first search for (typename MatrixType::InnerIterator it(A, jj); it; ++it) { - krow = it.row(); - dense_col(krow) = it.value(); - kmark = marker(krow); + int krow = it.row(); + dense_col(krow) = it.value(); + + int kmark = marker(krow); if (kmark == jj) continue; // krow visited before, go to the next nonzero - // For each unmarked krow of jj - marker(krow) = jj; - kperm = perm_r(krow); - if (kperm == IND_EMPTY ) { - // krow is in L : place it in structure of L(*, jj) - panel_lsub(nextl_col++) = krow; // krow is indexed into A - } - else - { - // krow is in U : if its supernode-representative krep - // has been explored, update repfnz(*) - krep = xsup(supno(kperm)+1) - 1; - myfnz = repfnz_col(krep); - - if (myfnz != IND_EMPTY ) - { - // Representative visited before - if (myfnz > kperm ) repfnz_col(krep) = kperm; - - } - else - { - // Otherwise, perform dfs starting at krep - oldrep = IND_EMPTY; - parent(krep) = oldrep; - repfnz_col(krep) = kperm; - xdfs = xlsub(krep); - maxdfs = xprune(krep); - - do - { - // For each unmarked kchild of krep - while (xdfs < maxdfs) - { - kchild = lsub(xdfs); - xdfs++; - chmark = marker(kchild); - - if (chmark != jj ) - { - marker(kchild) = jj; - chperm = perm_r(kchild); - - if (chperm == IND_EMPTY) - { - // case kchild is in L: place it in L(*, j) - panel_lsub(nextl_col++) = kchild; - } - else - { - // case kchild is in U : - // chrep = its supernode-rep. If its rep has been explored, - // update its repfnz(*) - chrep = xsup(supno(chperm)+1) - 1; - myfnz = repfnz_col(chrep); - - if (myfnz != IND_EMPTY) - { // Visited before - if (myfnz > chperm) - repfnz_col(chrep) = chperm; - } - else - { // Cont. dfs at snode-rep of kchild - xplore(krep) = xdfs; - oldrep = krep; - krep = chrep; // Go deeper down G(L) - parent(krep) = oldrep; - repfnz_col(krep) = chperm; - xdfs = xlsub(krep); - maxdfs = xprune(krep); - - } // end if myfnz != -1 - } // end if chperm == -1 - - } // end if chmark !=jj - } // end while xdfs < maxdfs - - // krow has no more unexplored nbrs : - // Place snode-rep krep in postorder DFS, if this - // segment is seen for the first time. (Note that - // "repfnz(krep)" may change later.) - // Baktrack dfs to its parent - if (marker1(krep) < jcol ) - { - segrep(nseg) = krep; - ++nseg; - marker1(krep) = jj; - } - - kpar = parent(krep); // Pop recursion, mimic recursion - if (kpar == IND_EMPTY) - break; // dfs done - krep = kpar; - xdfs = xplore(krep); - maxdfs = xprune(krep); - - } while (kpar != IND_EMPTY); // Do until empty stack - - } // end if (myfnz = -1) - - } // end if (kperm == -1) - + LU_dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent, + xplore, glu, nextl_col, krow, traits); }// end for nonzeros in column jj } // end for column jj