1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 StorageIndex; 41 panel_dfs_traits(Index jcol, StorageIndex* marker) 42 : m_jcol(jcol), m_marker(marker) 43 {} 44 bool update_segrep(Index krep, StorageIndex 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 StorageIndex* m_marker; 57 }; 58 59 60 template <typename Scalar, typename StorageIndex> 61 template <typename Traits> 62 void 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, 66 Index& nextl_col, Index krow, Traits& traits 67 ) 68 { 69 70 StorageIndex kmark = marker(krow); 71 72 // For each unmarked krow of jj 73 marker(krow) = jj; 74 StorageIndex 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++) = StorageIndex(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 StorageIndex krep = glu.xsup(glu.supno(kperm)+1) - 1; 87 // First nonzero element in the current column: 88 StorageIndex 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 StorageIndex oldrep = emptyIdxLU; 100 parent(krep) = oldrep; 101 repfnz_col(krep) = kperm; 102 StorageIndex xdfs = glu.xlsub(krep); 103 Index maxdfs = xprune(krep); 104 105 StorageIndex kpar; 106 do 107 { 108 // For each unmarked kchild of krep 109 while (xdfs < maxdfs) 110 { 111 StorageIndex kchild = glu.lsub(xdfs); 112 xdfs++; 113 StorageIndex chmark = marker(kchild); 114 115 if (chmark != jj ) 116 { 117 marker(kchild) = jj; 118 StorageIndex 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 StorageIndex 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 StorageIndex> 219 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) 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 (StorageIndex jj = StorageIndex(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 StorageIndex 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 259
