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]column_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_COLUMN_DFS_H
31#define SPARSELU_COLUMN_DFS_H
32
33template <typename Scalar, typename Index> class SparseLUImpl;
34namespace Eigen {
35
36namespace internal {
37
38template<typename IndexVector, typename ScalarVector>
39struct column_dfs_traits : no_assignment_operator
40{
41  typedef typename ScalarVector::Scalar Scalar;
42  typedef typename IndexVector::Scalar Index;
43  column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, Index>::GlobalLU_t& glu, SparseLUImpl<Scalar, Index>& luImpl)
44   : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
45 {}
46  bool update_segrep(Index /*krep*/, Index /*jj*/)
47  {
48    return true;
49  }
50  void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
51  {
52    if (nextl >= m_glu.nzlmax)
53      m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
54    if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU;
55  }
56  enum { ExpandMem = true };
57
58  Index m_jcol;
59  Index& m_jsuper_ref;
60  typename SparseLUImpl<Scalar, Index>::GlobalLU_t& m_glu;
61  SparseLUImpl<Scalar, Index>& m_luImpl;
62};
63
64
65/**
66 * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
67 *
68 * A supernode representative is the last column of a supernode.
69 * The nonzeros in U[*,j] are segments that end at supernodes representatives.
70 * The routine returns a list of the supernodal representatives
71 * in topological order of the dfs that generates them.
72 * The location of the first nonzero in each supernodal segment
73 * (supernodal entry location) is also returned.
74 *
75 * \param m number of rows in the matrix
76 * \param jcol Current column
77 * \param perm_r Row permutation
78 * \param maxsuper  Maximum number of column allowed in a supernode
79 * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
80 * \param lsub_col defines the rhs vector to start the dfs
81 * \param [in,out] segrep Segment representatives - new segments appended
82 * \param repfnz  First nonzero location in each row
83 * \param xprune
84 * \param marker  marker[i] == jj, if i was visited during dfs of current column jj;
85 * \param parent
86 * \param xplore working array
87 * \param glu global LU data
88 * \return 0 success
89 *         > 0 number of bytes allocated when run out of space
90 *
91 */
92template <typename Scalar, typename Index>
93Index SparseLUImpl<Scalar,Index>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,  BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
94{
95
96  Index jsuper = glu.supno(jcol);
97  Index nextl = glu.xlsub(jcol);
98  VectorBlock<IndexVector> marker2(marker, 2*m, m);
99
100
101  column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
102
103  // For each nonzero in A(*,jcol) do dfs
104  for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
105  {
106    Index krow = lsub_col(k);
107    lsub_col(k) = emptyIdxLU;
108    Index kmark = marker2(krow);
109
110    // krow was visited before, go to the next nonz;
111    if (kmark == jcol) continue;
112
113    dfs_kernel(jcol, perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
114                   xplore, glu, nextl, krow, traits);
115  } // for each nonzero ...
116
117  Index fsupc, jptr, jm1ptr, ito, ifrom, istop;
118  Index nsuper = glu.supno(jcol);
119  Index jcolp1 = jcol + 1;
120  Index jcolm1 = jcol - 1;
121
122  // check to see if j belongs in the same supernode as j-1
123  if ( jcol == 0 )
124  { // Do nothing for column 0
125    nsuper = glu.supno(0) = 0 ;
126  }
127  else
128  {
129    fsupc = glu.xsup(nsuper);
130    jptr = glu.xlsub(jcol); // Not yet compressed
131    jm1ptr = glu.xlsub(jcolm1);
132
133    // Use supernodes of type T2 : see SuperLU paper
134    if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
135
136    // Make sure the number of columns in a supernode doesn't
137    // exceed threshold
138    if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU;
139
140    /* If jcol starts a new supernode, reclaim storage space in
141     * glu.lsub from previous supernode. Note we only store
142     * the subscript set of the first and last columns of
143     * a supernode. (first for num values, last for pruning)
144     */
145    if (jsuper == emptyIdxLU)
146    { // starts a new supernode
147      if ( (fsupc < jcolm1-1) )
148      { // >= 3 columns in nsuper
149        ito = glu.xlsub(fsupc+1);
150        glu.xlsub(jcolm1) = ito;
151        istop = ito + jptr - jm1ptr;
152        xprune(jcolm1) = istop; // intialize xprune(jcol-1)
153        glu.xlsub(jcol) = istop;
154
155        for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
156          glu.lsub(ito) = glu.lsub(ifrom);
157        nextl = ito;  // = istop + length(jcol)
158      }
159      nsuper++;
160      glu.supno(jcol) = nsuper;
161    } // if a new supernode
162  } // end else:  jcol > 0
163
164  // Tidy up the pointers before exit
165  glu.xsup(nsuper+1) = jcolp1;
166  glu.supno(jcolp1) = nsuper;
167  xprune(jcol) = nextl;  // Intialize upper bound for pruning
168  glu.xlsub(jcolp1) = nextl;
169
170  return 0;
171}
172
173} // end namespace internal
174
175} // end namespace Eigen
176
177#endif
178