SparseLU_panel_dfs.h
Go to the documentation of this file.
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>
39 {
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 };
56  StorageIndex* m_marker;
57 };
58 
59 
60 template <typename Scalar, typename StorageIndex>
61 template <typename Traits>
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 
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
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:33
panel_dfs_traits(Index jcol, StorageIndex *marker)
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
void 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)
Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
void mem_expand(IndexVector &, Index, Index)
bool update_segrep(Index krep, StorageIndex jj)
Expression of a fixed-size or dynamic-size sub-vector.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
RowVector3d w
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:192
void 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)


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:25