| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| /* |
| |
| * NOTE: This file is the modified version of [s,d,c,z]panel_dfs.c file in SuperLU |
| |
| * -- SuperLU routine (version 2.0) -- |
| * Univ. of California Berkeley, Xerox Palo Alto Research Center, |
| * and Lawrence Berkeley National Lab. |
| * November 15, 1997 |
| * |
| * Copyright (c) 1994 by Xerox Corporation. All rights reserved. |
| * |
| * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
| * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
| * |
| * Permission is hereby granted to use or copy this program for any |
| * purpose, provided the above notices are retained on all copies. |
| * Permission to modify the code and to distribute modified code is |
| * granted, provided the above notices are retained, and a notice that |
| * the code was modified is included with the above copyright notice. |
| */ |
| #ifndef SPARSELU_PANEL_DFS_H |
| #define SPARSELU_PANEL_DFS_H |
| |
| // IWYU pragma: private |
| #include "./InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| namespace internal { |
| |
| template <typename IndexVector> |
| struct panel_dfs_traits { |
| typedef typename IndexVector::Scalar StorageIndex; |
| panel_dfs_traits(Index jcol, StorageIndex* marker) : m_jcol(jcol), m_marker(marker) {} |
| bool update_segrep(Index krep, StorageIndex jj) { |
| if (m_marker[krep] < m_jcol) { |
| m_marker[krep] = jj; |
| return true; |
| } |
| return false; |
| } |
| void mem_expand(IndexVector& /*glu.lsub*/, Index /*nextl*/, Index /*chmark*/) {} |
| enum { ExpandMem = false }; |
| Index m_jcol; |
| StorageIndex* m_marker; |
| }; |
| |
| template <typename Scalar, typename StorageIndex> |
| template <typename Traits> |
| void SparseLUImpl<Scalar, StorageIndex>::dfs_kernel(const StorageIndex jj, IndexVector& perm_r, Index& nseg, |
| IndexVector& panel_lsub, IndexVector& segrep, |
| Ref<IndexVector> repfnz_col, IndexVector& xprune, |
| Ref<IndexVector> marker, IndexVector& parent, IndexVector& xplore, |
| GlobalLU_t& glu, Index& nextl_col, Index krow, Traits& traits) { |
| StorageIndex kmark = marker(krow); |
| |
| // For each unmarked krow of jj |
| marker(krow) = jj; |
| StorageIndex kperm = perm_r(krow); |
| if (kperm == emptyIdxLU) { |
| // krow is in L : place it in structure of L(*, jj) |
| panel_lsub(nextl_col++) = StorageIndex(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 |
| StorageIndex krep = glu.xsup(glu.supno(kperm) + 1) - 1; |
| // First nonzero element in the current column: |
| StorageIndex myfnz = repfnz_col(krep); |
| |
| if (myfnz != emptyIdxLU) { |
| // Representative visited before |
| if (myfnz > kperm) repfnz_col(krep) = kperm; |
| |
| } else { |
| // Otherwise, perform dfs starting at krep |
| StorageIndex oldrep = emptyIdxLU; |
| parent(krep) = oldrep; |
| repfnz_col(krep) = kperm; |
| StorageIndex xdfs = glu.xlsub(krep); |
| Index maxdfs = xprune(krep); |
| |
| StorageIndex kpar; |
| do { |
| // For each unmarked kchild of krep |
| while (xdfs < maxdfs) { |
| StorageIndex kchild = glu.lsub(xdfs); |
| xdfs++; |
| StorageIndex chmark = marker(kchild); |
| |
| if (chmark != jj) { |
| marker(kchild) = jj; |
| StorageIndex chperm = perm_r(kchild); |
| |
| if (chperm == emptyIdxLU) { |
| // 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(*) |
| StorageIndex chrep = glu.xsup(glu.supno(chperm) + 1) - 1; |
| myfnz = repfnz_col(chrep); |
| |
| if (myfnz != emptyIdxLU) { // 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 = glu.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 == emptyIdxLU) break; // dfs done |
| krep = kpar; |
| xdfs = xplore(krep); |
| maxdfs = xprune(krep); |
| |
| } while (kpar != emptyIdxLU); // 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) |
| * |
| * A supernode representative is the last column of a supernode. |
| * The nonzeros in U[*,j] are segments that end at supernodes representatives |
| * |
| * The routine returns a list of the supernodal representatives |
| * in topological order of the dfs that generates them. This list is |
| * a superset of the topological order of each individual column within |
| * the panel. |
| * The location of the first nonzero in each supernodal segment |
| * (supernodal entry location) is also returned. Each column has |
| * a separate list for this purpose. |
| * |
| * Two markers arrays are used for dfs : |
| * marker[i] == jj, if i was visited during dfs of current column jj; |
| * marker1[i] >= jcol, if i was visited by earlier columns in this panel; |
| * |
| * \param[in] m number of rows in the matrix |
| * \param[in] w Panel size |
| * \param[in] jcol Starting column of the panel |
| * \param[in] A Input matrix in column-major storage |
| * \param[in] perm_r Row permutation |
| * \param[out] nseg Number of U segments |
| * \param[out] dense Accumulate the column vectors of the panel |
| * \param[out] panel_lsub Subscripts of the row in the panel |
| * \param[out] segrep Segment representative i.e first nonzero row of each segment |
| * \param[out] repfnz First nonzero location in each row |
| * \param[out] xprune The pruned elimination tree |
| * \param[out] marker work vector |
| * \param parent The elimination tree |
| * \param xplore work vector |
| * \param glu The global data structure |
| * |
| */ |
| |
| template <typename Scalar, typename StorageIndex> |
| void 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) { |
| Index nextl_col; // Next available position in panel_lsub[*,jj] |
| |
| // Initialize pointers |
| VectorBlock<IndexVector> marker1(marker, m, m); |
| nseg = 0; |
| |
| panel_dfs_traits<IndexVector> traits(jcol, marker1.data()); |
| |
| // For each column in the panel |
| for (StorageIndex jj = StorageIndex(jcol); jj < jcol + w; jj++) { |
| nextl_col = (jj - jcol) * m; |
| |
| VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row |
| VectorBlock<ScalarVector> dense_col(dense, nextl_col, m); // Accumulate a column vector here |
| |
| // For each nnz in A[*, jj] do depth first search |
| for (typename MatrixType::InnerIterator it(A, jj); it; ++it) { |
| Index krow = it.row(); |
| dense_col(krow) = it.value(); |
| |
| StorageIndex kmark = marker(krow); |
| if (kmark == jj) continue; // krow visited before, go to the next nonzero |
| |
| 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 |
| } |
| |
| } // end namespace internal |
| } // end namespace Eigen |
| |
| #endif // SPARSELU_PANEL_DFS_H |
