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
13 * NOTE: This file is the modified version of sp_coletree.c file in SuperLU
14
15 * -- SuperLU routine (version 3.1) --
16 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
17 * and Lawrence Berkeley National Lab.
18 * August 1, 2008
19 *
20 * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
21 *
22 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
23 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
24 *
25 * Permission is hereby granted to use or copy this program for any
26 * purpose, provided the above notices are retained on all copies.
27 * Permission to modify the code and to distribute modified code is
28 * granted, provided the above notices are retained, and a notice that
29 * the code was modified is included with the above copyright notice.
30 */
31 #ifndef SPARSE_COLETREE_H
32 #define SPARSE_COLETREE_H
33
34 namespace Eigen {
35
36 namespace internal {
37
38 /** Find the root of the tree/set containing the vertex i : Use Path halving */
39 template<typename Index, typename IndexVector>
etree_find(Index i,IndexVector & pp)40 Index etree_find (Index i, IndexVector& pp)
41 {
42 Index p = pp(i); // Parent
43 Index gp = pp(p); // Grand parent
44 while (gp != p)
45 {
46 pp(i) = gp; // Parent pointer on find path is changed to former grand parent
47 i = gp;
48 p = pp(i);
49 gp = pp(p);
50 }
51 return p;
52 }
53
54 /** Compute the column elimination tree of a sparse matrix
55 * \param mat The matrix in column-major format.
56 * \param parent The elimination tree
57 * \param firstRowElt The column index of the first element in each row
58 * \param perm The permutation to apply to the column of \b mat
59 */
60 template <typename MatrixType, typename IndexVector>
61 int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::StorageIndex *perm=0)
62 {
63 typedef typename MatrixType::StorageIndex StorageIndex;
64 StorageIndex nc = convert_index<StorageIndex>(mat.cols()); // Number of columns
65 StorageIndex m = convert_index<StorageIndex>(mat.rows());
66 StorageIndex diagSize = (std::min)(nc,m);
67 IndexVector root(nc); // root of subtree of etree
68 root.setZero();
69 IndexVector pp(nc); // disjoint sets
70 pp.setZero(); // Initialize disjoint sets
71 parent.resize(mat.cols());
72 //Compute first nonzero column in each row
73 firstRowElt.resize(m);
74 firstRowElt.setConstant(nc);
75 firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
76 bool found_diag;
77 for (StorageIndex col = 0; col < nc; col++)
78 {
79 StorageIndex pcol = col;
80 if(perm) pcol = perm[col];
81 for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
82 {
83 Index row = it.row();
84 firstRowElt(row) = (std::min)(firstRowElt(row), col);
85 }
86 }
87 /* Compute etree by Liu's algorithm for symmetric matrices,
88 except use (firstRowElt[r],c) in place of an edge (r,c) of A.
89 Thus each row clique in A'*A is replaced by a star
90 centered at its first vertex, which has the same fill. */
91 StorageIndex rset, cset, rroot;
92 for (StorageIndex col = 0; col < nc; col++)
93 {
94 found_diag = col>=m;
95 pp(col) = col;
96 cset = col;
97 root(cset) = col;
98 parent(col) = nc;
99 /* The diagonal element is treated here even if it does not exist in the matrix
100 * hence the loop is executed once more */
101 StorageIndex pcol = col;
102 if(perm) pcol = perm[col];
103 for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
104 { // A sequence of interleaved find and union is performed
105 Index i = col;
106 if(it) i = it.index();
107 if (i == col) found_diag = true;
108
109 StorageIndex row = firstRowElt(i);
110 if (row >= col) continue;
111 rset = internal::etree_find(row, pp); // Find the name of the set containing row
112 rroot = root(rset);
113 if (rroot != col)
114 {
115 parent(rroot) = col;
116 pp(cset) = rset;
117 cset = rset;
118 root(cset) = col;
119 }
120 }
121 }
122 return 0;
123 }
124
125 /**
126 * Depth-first search from vertex n. No recursion.
127 * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
128 */
129 template <typename IndexVector>
nr_etdfs(typename IndexVector::Scalar n,IndexVector & parent,IndexVector & first_kid,IndexVector & next_kid,IndexVector & post,typename IndexVector::Scalar postnum)130 void nr_etdfs (typename IndexVector::Scalar n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, typename IndexVector::Scalar postnum)
131 {
132 typedef typename IndexVector::Scalar StorageIndex;
133 StorageIndex current = n, first, next;
134 while (postnum != n)
135 {
136 // No kid for the current node
137 first = first_kid(current);
138
139 // no kid for the current node
140 if (first == -1)
141 {
142 // Numbering this node because it has no kid
143 post(current) = postnum++;
144
145 // looking for the next kid
146 next = next_kid(current);
147 while (next == -1)
148 {
149 // No more kids : back to the parent node
150 current = parent(current);
151 // numbering the parent node
152 post(current) = postnum++;
153
154 // Get the next kid
155 next = next_kid(current);
156 }
157 // stopping criterion
158 if (postnum == n+1) return;
159
160 // Updating current node
161 current = next;
162 }
163 else
164 {
165 current = first;
166 }
167 }
168 }
169
170
171 /**
172 * \brief Post order a tree
173 * \param n the number of nodes
174 * \param parent Input tree
175 * \param post postordered tree
176 */
177 template <typename IndexVector>
treePostorder(typename IndexVector::Scalar n,IndexVector & parent,IndexVector & post)178 void treePostorder(typename IndexVector::Scalar n, IndexVector& parent, IndexVector& post)
179 {
180 typedef typename IndexVector::Scalar StorageIndex;
181 IndexVector first_kid, next_kid; // Linked list of children
182 StorageIndex postnum;
183 // Allocate storage for working arrays and results
184 first_kid.resize(n+1);
185 next_kid.setZero(n+1);
186 post.setZero(n+1);
187
188 // Set up structure describing children
189 first_kid.setConstant(-1);
190 for (StorageIndex v = n-1; v >= 0; v--)
191 {
192 StorageIndex dad = parent(v);
193 next_kid(v) = first_kid(dad);
194 first_kid(dad) = v;
195 }
196
197 // Depth-first search from dummy root vertex #n
198 postnum = 0;
199 internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
200 }
201
202 } // end namespace internal
203
204 } // end namespace Eigen
205
206 #endif // SPARSE_COLETREE_H
207