| 1 | // This file is part of Eigen, a lightweight C++ template library
|
|---|
| 2 | // for linear algebra.
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
|
|---|
| 5 | //
|
|---|
| 6 | // This Source Code Form is subject to the terms of the Mozilla
|
|---|
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
|---|
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|---|
| 9 |
|
|---|
| 10 | /*
|
|---|
| 11 |
|
|---|
| 12 | * NOTE: This file is the modified version of [s,d,c,z]panel_dfs.c file in SuperLU
|
|---|
| 13 |
|
|---|
| 14 | * -- SuperLU routine (version 2.0) --
|
|---|
| 15 | * Univ. of California Berkeley, Xerox Palo Alto Research Center,
|
|---|
| 16 | * and Lawrence Berkeley National Lab.
|
|---|
| 17 | * November 15, 1997
|
|---|
| 18 | *
|
|---|
| 19 | * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
|
|---|
| 20 | *
|
|---|
| 21 | * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
|
|---|
| 22 | * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
|
|---|
| 23 | *
|
|---|
| 24 | * Permission is hereby granted to use or copy this program for any
|
|---|
| 25 | * purpose, provided the above notices are retained on all copies.
|
|---|
| 26 | * Permission to modify the code and to distribute modified code is
|
|---|
| 27 | * granted, provided the above notices are retained, and a notice that
|
|---|
| 28 | * the code was modified is included with the above copyright notice.
|
|---|
| 29 | */
|
|---|
| 30 | #ifndef SPARSELU_PANEL_DFS_H
|
|---|
| 31 | #define SPARSELU_PANEL_DFS_H
|
|---|
| 32 |
|
|---|
| 33 | namespace Eigen {
|
|---|
| 34 |
|
|---|
| 35 | namespace internal {
|
|---|
| 36 |
|
|---|
| 37 | template<typename IndexVector>
|
|---|
| 38 | struct panel_dfs_traits
|
|---|
| 39 | {
|
|---|
| 40 | typedef typename IndexVector::Scalar Index;
|
|---|
| 41 | panel_dfs_traits(Index jcol, Index* marker)
|
|---|
| 42 | : m_jcol(jcol), m_marker(marker)
|
|---|
| 43 | {}
|
|---|
| 44 | bool update_segrep(Index krep, Index jj)
|
|---|
| 45 | {
|
|---|
| 46 | if(m_marker[krep]<m_jcol)
|
|---|
| 47 | {
|
|---|
| 48 | m_marker[krep] = jj;
|
|---|
| 49 | return true;
|
|---|
| 50 | }
|
|---|
| 51 | return false;
|
|---|
| 52 | }
|
|---|
| 53 | void mem_expand(IndexVector& /*glu.lsub*/, Index /*nextl*/, Index /*chmark*/) {}
|
|---|
| 54 | enum { ExpandMem = false };
|
|---|
| 55 | Index m_jcol;
|
|---|
| 56 | Index* m_marker;
|
|---|
| 57 | };
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 | template <typename Scalar, typename Index>
|
|---|
| 61 | template <typename Traits>
|
|---|
| 62 | void SparseLUImpl<Scalar,Index>::dfs_kernel(const Index 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,
|
|---|
| 66 | Index& nextl_col, Index krow, Traits& traits
|
|---|
| 67 | )
|
|---|
| 68 | {
|
|---|
| 69 |
|
|---|
| 70 | Index kmark = marker(krow);
|
|---|
| 71 |
|
|---|
| 72 | // For each unmarked krow of jj
|
|---|
| 73 | marker(krow) = jj;
|
|---|
| 74 | Index kperm = perm_r(krow);
|
|---|
| 75 | if (kperm == emptyIdxLU ) {
|
|---|
| 76 | // krow is in L : place it in structure of L(*, jj)
|
|---|
| 77 | panel_lsub(nextl_col++) = krow; // krow is indexed into A
|
|---|
| 78 |
|
|---|
| 79 | traits.mem_expand(panel_lsub, nextl_col, kmark);
|
|---|
| 80 | }
|
|---|
| 81 | else
|
|---|
| 82 | {
|
|---|
| 83 | // krow is in U : if its supernode-representative krep
|
|---|
| 84 | // has been explored, update repfnz(*)
|
|---|
| 85 | // krep = supernode representative of the current row
|
|---|
| 86 | Index krep = glu.xsup(glu.supno(kperm)+1) - 1;
|
|---|
| 87 | // First nonzero element in the current column:
|
|---|
| 88 | Index myfnz = repfnz_col(krep);
|
|---|
| 89 |
|
|---|
| 90 | if (myfnz != emptyIdxLU )
|
|---|
| 91 | {
|
|---|
| 92 | // Representative visited before
|
|---|
| 93 | if (myfnz > kperm ) repfnz_col(krep) = kperm;
|
|---|
| 94 |
|
|---|
| 95 | }
|
|---|
| 96 | else
|
|---|
| 97 | {
|
|---|
| 98 | // Otherwise, perform dfs starting at krep
|
|---|
| 99 | Index oldrep = emptyIdxLU;
|
|---|
| 100 | parent(krep) = oldrep;
|
|---|
| 101 | repfnz_col(krep) = kperm;
|
|---|
| 102 | Index xdfs = glu.xlsub(krep);
|
|---|
| 103 | Index maxdfs = xprune(krep);
|
|---|
| 104 |
|
|---|
| 105 | Index kpar;
|
|---|
| 106 | do
|
|---|
| 107 | {
|
|---|
| 108 | // For each unmarked kchild of krep
|
|---|
| 109 | while (xdfs < maxdfs)
|
|---|
| 110 | {
|
|---|
| 111 | Index kchild = glu.lsub(xdfs);
|
|---|
| 112 | xdfs++;
|
|---|
| 113 | Index chmark = marker(kchild);
|
|---|
| 114 |
|
|---|
| 115 | if (chmark != jj )
|
|---|
| 116 | {
|
|---|
| 117 | marker(kchild) = jj;
|
|---|
| 118 | Index chperm = perm_r(kchild);
|
|---|
| 119 |
|
|---|
| 120 | if (chperm == emptyIdxLU)
|
|---|
| 121 | {
|
|---|
| 122 | // case kchild is in L: place it in L(*, j)
|
|---|
| 123 | panel_lsub(nextl_col++) = kchild;
|
|---|
| 124 | traits.mem_expand(panel_lsub, nextl_col, chmark);
|
|---|
| 125 | }
|
|---|
| 126 | else
|
|---|
| 127 | {
|
|---|
| 128 | // case kchild is in U :
|
|---|
| 129 | // chrep = its supernode-rep. If its rep has been explored,
|
|---|
| 130 | // update its repfnz(*)
|
|---|
| 131 | Index chrep = glu.xsup(glu.supno(chperm)+1) - 1;
|
|---|
| 132 | myfnz = repfnz_col(chrep);
|
|---|
| 133 |
|
|---|
| 134 | if (myfnz != emptyIdxLU)
|
|---|
| 135 | { // Visited before
|
|---|
| 136 | if (myfnz > chperm)
|
|---|
| 137 | repfnz_col(chrep) = chperm;
|
|---|
| 138 | }
|
|---|
| 139 | else
|
|---|
| 140 | { // Cont. dfs at snode-rep of kchild
|
|---|
| 141 | xplore(krep) = xdfs;
|
|---|
| 142 | oldrep = krep;
|
|---|
| 143 | krep = chrep; // Go deeper down G(L)
|
|---|
| 144 | parent(krep) = oldrep;
|
|---|
| 145 | repfnz_col(krep) = chperm;
|
|---|
| 146 | xdfs = glu.xlsub(krep);
|
|---|
| 147 | maxdfs = xprune(krep);
|
|---|
| 148 |
|
|---|
| 149 | } // end if myfnz != -1
|
|---|
| 150 | } // end if chperm == -1
|
|---|
| 151 |
|
|---|
| 152 | } // end if chmark !=jj
|
|---|
| 153 | } // end while xdfs < maxdfs
|
|---|
| 154 |
|
|---|
| 155 | // krow has no more unexplored nbrs :
|
|---|
| 156 | // Place snode-rep krep in postorder DFS, if this
|
|---|
| 157 | // segment is seen for the first time. (Note that
|
|---|
| 158 | // "repfnz(krep)" may change later.)
|
|---|
| 159 | // Baktrack dfs to its parent
|
|---|
| 160 | if(traits.update_segrep(krep,jj))
|
|---|
| 161 | //if (marker1(krep) < jcol )
|
|---|
| 162 | {
|
|---|
| 163 | segrep(nseg) = krep;
|
|---|
| 164 | ++nseg;
|
|---|
| 165 | //marker1(krep) = jj;
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | kpar = parent(krep); // Pop recursion, mimic recursion
|
|---|
| 169 | if (kpar == emptyIdxLU)
|
|---|
| 170 | break; // dfs done
|
|---|
| 171 | krep = kpar;
|
|---|
| 172 | xdfs = xplore(krep);
|
|---|
| 173 | maxdfs = xprune(krep);
|
|---|
| 174 |
|
|---|
| 175 | } while (kpar != emptyIdxLU); // Do until empty stack
|
|---|
| 176 |
|
|---|
| 177 | } // end if (myfnz = -1)
|
|---|
| 178 |
|
|---|
| 179 | } // end if (kperm == -1)
|
|---|
| 180 | }
|
|---|
| 181 |
|
|---|
| 182 | /**
|
|---|
| 183 | * \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
|
|---|
| 184 | *
|
|---|
| 185 | * A supernode representative is the last column of a supernode.
|
|---|
| 186 | * The nonzeros in U[*,j] are segments that end at supernodes representatives
|
|---|
| 187 | *
|
|---|
| 188 | * The routine returns a list of the supernodal representatives
|
|---|
| 189 | * in topological order of the dfs that generates them. This list is
|
|---|
| 190 | * a superset of the topological order of each individual column within
|
|---|
| 191 | * the panel.
|
|---|
| 192 | * The location of the first nonzero in each supernodal segment
|
|---|
| 193 | * (supernodal entry location) is also returned. Each column has
|
|---|
| 194 | * a separate list for this purpose.
|
|---|
| 195 | *
|
|---|
| 196 | * Two markers arrays are used for dfs :
|
|---|
| 197 | * marker[i] == jj, if i was visited during dfs of current column jj;
|
|---|
| 198 | * marker1[i] >= jcol, if i was visited by earlier columns in this panel;
|
|---|
| 199 | *
|
|---|
| 200 | * \param[in] m number of rows in the matrix
|
|---|
| 201 | * \param[in] w Panel size
|
|---|
| 202 | * \param[in] jcol Starting column of the panel
|
|---|
| 203 | * \param[in] A Input matrix in column-major storage
|
|---|
| 204 | * \param[in] perm_r Row permutation
|
|---|
| 205 | * \param[out] nseg Number of U segments
|
|---|
| 206 | * \param[out] dense Accumulate the column vectors of the panel
|
|---|
| 207 | * \param[out] panel_lsub Subscripts of the row in the panel
|
|---|
| 208 | * \param[out] segrep Segment representative i.e first nonzero row of each segment
|
|---|
| 209 | * \param[out] repfnz First nonzero location in each row
|
|---|
| 210 | * \param[out] xprune The pruned elimination tree
|
|---|
| 211 | * \param[out] marker work vector
|
|---|
| 212 | * \param parent The elimination tree
|
|---|
| 213 | * \param xplore work vector
|
|---|
| 214 | * \param glu The global data structure
|
|---|
| 215 | *
|
|---|
| 216 | */
|
|---|
| 217 |
|
|---|
| 218 | template <typename Scalar, typename Index>
|
|---|
| 219 | void SparseLUImpl<Scalar,Index>::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)
|
|---|
| 220 | {
|
|---|
| 221 | Index nextl_col; // Next available position in panel_lsub[*,jj]
|
|---|
| 222 |
|
|---|
| 223 | // Initialize pointers
|
|---|
| 224 | VectorBlock<IndexVector> marker1(marker, m, m);
|
|---|
| 225 | nseg = 0;
|
|---|
| 226 |
|
|---|
| 227 | panel_dfs_traits<IndexVector> traits(jcol, marker1.data());
|
|---|
| 228 |
|
|---|
| 229 | // For each column in the panel
|
|---|
| 230 | for (Index jj = jcol; jj < jcol + w; jj++)
|
|---|
| 231 | {
|
|---|
| 232 | nextl_col = (jj - jcol) * m;
|
|---|
| 233 |
|
|---|
| 234 | VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row
|
|---|
| 235 | VectorBlock<ScalarVector> dense_col(dense,nextl_col, m); // Accumulate a column vector here
|
|---|
| 236 |
|
|---|
| 237 |
|
|---|
| 238 | // For each nnz in A[*, jj] do depth first search
|
|---|
| 239 | for (typename MatrixType::InnerIterator it(A, jj); it; ++it)
|
|---|
| 240 | {
|
|---|
| 241 | Index krow = it.row();
|
|---|
| 242 | dense_col(krow) = it.value();
|
|---|
| 243 |
|
|---|
| 244 | Index kmark = marker(krow);
|
|---|
| 245 | if (kmark == jj)
|
|---|
| 246 | continue; // krow visited before, go to the next nonzero
|
|---|
| 247 |
|
|---|
| 248 | dfs_kernel(jj, perm_r, nseg, panel_lsub, segrep, repfnz_col, xprune, marker, parent,
|
|---|
| 249 | xplore, glu, nextl_col, krow, traits);
|
|---|
| 250 | }// end for nonzeros in column jj
|
|---|
| 251 |
|
|---|
| 252 | } // end for column jj
|
|---|
| 253 | }
|
|---|
| 254 |
|
|---|
| 255 | } // end namespace internal
|
|---|
| 256 | } // end namespace Eigen
|
|---|
| 257 |
|
|---|
| 258 | #endif // SPARSELU_PANEL_DFS_H
|
|---|
