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 
33 template <typename Scalar, typename StorageIndex> class SparseLUImpl;
34 namespace Eigen {
35 
36 namespace internal {
37 
38 template<typename IndexVector, typename ScalarVector>
39 struct column_dfs_traits : no_assignment_operator
40 {
41   typedef typename ScalarVector::Scalar Scalar;
42   typedef typename IndexVector::Scalar StorageIndex;
column_dfs_traitscolumn_dfs_traits43   column_dfs_traits(Index jcol, Index& jsuper, typename SparseLUImpl<Scalar, StorageIndex>::GlobalLU_t& glu, SparseLUImpl<Scalar, StorageIndex>& luImpl)
44    : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
45  {}
update_segrepcolumn_dfs_traits46   bool update_segrep(Index /*krep*/, Index /*jj*/)
47   {
48     return true;
49   }
mem_expandcolumn_dfs_traits50   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, StorageIndex>::GlobalLU_t& m_glu;
61   SparseLUImpl<Scalar, StorageIndex>& 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  */
92 template <typename Scalar, typename StorageIndex>
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)93 Index SparseLUImpl<Scalar,StorageIndex>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg,
94                                                     BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune,
95                                                     IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
96 {
97 
98   Index jsuper = glu.supno(jcol);
99   Index nextl = glu.xlsub(jcol);
100   VectorBlock<IndexVector> marker2(marker, 2*m, m);
101 
102 
103   column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
104 
105   // For each nonzero in A(*,jcol) do dfs
106   for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
107   {
108     Index krow = lsub_col(k);
109     lsub_col(k) = emptyIdxLU;
110     Index kmark = marker2(krow);
111 
112     // krow was visited before, go to the next nonz;
113     if (kmark == jcol) continue;
114 
115     dfs_kernel(StorageIndex(jcol), perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
116                    xplore, glu, nextl, krow, traits);
117   } // for each nonzero ...
118 
119   Index fsupc;
120   StorageIndex nsuper = glu.supno(jcol);
121   StorageIndex jcolp1 = StorageIndex(jcol) + 1;
122   Index jcolm1 = jcol - 1;
123 
124   // check to see if j belongs in the same supernode as j-1
125   if ( jcol == 0 )
126   { // Do nothing for column 0
127     nsuper = glu.supno(0) = 0 ;
128   }
129   else
130   {
131     fsupc = glu.xsup(nsuper);
132     StorageIndex jptr = glu.xlsub(jcol); // Not yet compressed
133     StorageIndex jm1ptr = glu.xlsub(jcolm1);
134 
135     // Use supernodes of type T2 : see SuperLU paper
136     if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
137 
138     // Make sure the number of columns in a supernode doesn't
139     // exceed threshold
140     if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU;
141 
142     /* If jcol starts a new supernode, reclaim storage space in
143      * glu.lsub from previous supernode. Note we only store
144      * the subscript set of the first and last columns of
145      * a supernode. (first for num values, last for pruning)
146      */
147     if (jsuper == emptyIdxLU)
148     { // starts a new supernode
149       if ( (fsupc < jcolm1-1) )
150       { // >= 3 columns in nsuper
151         StorageIndex ito = glu.xlsub(fsupc+1);
152         glu.xlsub(jcolm1) = ito;
153         StorageIndex istop = ito + jptr - jm1ptr;
154         xprune(jcolm1) = istop; // intialize xprune(jcol-1)
155         glu.xlsub(jcol) = istop;
156 
157         for (StorageIndex ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
158           glu.lsub(ito) = glu.lsub(ifrom);
159         nextl = ito;  // = istop + length(jcol)
160       }
161       nsuper++;
162       glu.supno(jcol) = nsuper;
163     } // if a new supernode
164   } // end else:  jcol > 0
165 
166   // Tidy up the pointers before exit
167   glu.xsup(nsuper+1) = jcolp1;
168   glu.supno(jcolp1) = nsuper;
169   xprune(jcol) = StorageIndex(nextl);  // Intialize upper bound for pruning
170   glu.xlsub(jcolp1) = StorageIndex(nextl);
171 
172   return 0;
173 }
174 
175 } // end namespace internal
176 
177 } // end namespace Eigen
178 
179 #endif
180