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
|
---|