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
33namespace Eigen {
34
35namespace internal {
36
37template<typename IndexVector>
38struct 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
60template <typename Scalar, typename Index>
61template <typename Traits>
62void 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
218template <typename Scalar, typename Index>
219void 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
259