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>
|
---|
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::Index *perm=0)
|
---|
62 | {
|
---|
63 | typedef typename MatrixType::Index Index;
|
---|
64 | Index nc = mat.cols(); // Number of columns
|
---|
65 | Index m = mat.rows();
|
---|
66 | Index 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 | Index row,col;
|
---|
74 | firstRowElt.resize(m);
|
---|
75 | firstRowElt.setConstant(nc);
|
---|
76 | firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
|
---|
77 | bool found_diag;
|
---|
78 | for (col = 0; col < nc; col++)
|
---|
79 | {
|
---|
80 | Index pcol = col;
|
---|
81 | if(perm) pcol = perm[col];
|
---|
82 | for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
|
---|
83 | {
|
---|
84 | row = it.row();
|
---|
85 | firstRowElt(row) = (std::min)(firstRowElt(row), col);
|
---|
86 | }
|
---|
87 | }
|
---|
88 | /* Compute etree by Liu's algorithm for symmetric matrices,
|
---|
89 | except use (firstRowElt[r],c) in place of an edge (r,c) of A.
|
---|
90 | Thus each row clique in A'*A is replaced by a star
|
---|
91 | centered at its first vertex, which has the same fill. */
|
---|
92 | Index rset, cset, rroot;
|
---|
93 | for (col = 0; col < nc; col++)
|
---|
94 | {
|
---|
95 | found_diag = col>=m;
|
---|
96 | pp(col) = col;
|
---|
97 | cset = col;
|
---|
98 | root(cset) = col;
|
---|
99 | parent(col) = nc;
|
---|
100 | /* The diagonal element is treated here even if it does not exist in the matrix
|
---|
101 | * hence the loop is executed once more */
|
---|
102 | Index pcol = col;
|
---|
103 | if(perm) pcol = perm[col];
|
---|
104 | for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
|
---|
105 | { // A sequence of interleaved find and union is performed
|
---|
106 | Index i = col;
|
---|
107 | if(it) i = it.index();
|
---|
108 | if (i == col) found_diag = true;
|
---|
109 |
|
---|
110 | row = firstRowElt(i);
|
---|
111 | if (row >= col) continue;
|
---|
112 | rset = internal::etree_find(row, pp); // Find the name of the set containing row
|
---|
113 | rroot = root(rset);
|
---|
114 | if (rroot != col)
|
---|
115 | {
|
---|
116 | parent(rroot) = col;
|
---|
117 | pp(cset) = rset;
|
---|
118 | cset = rset;
|
---|
119 | root(cset) = col;
|
---|
120 | }
|
---|
121 | }
|
---|
122 | }
|
---|
123 | return 0;
|
---|
124 | }
|
---|
125 |
|
---|
126 | /**
|
---|
127 | * Depth-first search from vertex n. No recursion.
|
---|
128 | * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
|
---|
129 | */
|
---|
130 | template <typename Index, typename IndexVector>
|
---|
131 | void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum)
|
---|
132 | {
|
---|
133 | Index 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 Index, typename IndexVector>
|
---|
178 | void treePostorder(Index n, IndexVector& parent, IndexVector& post)
|
---|
179 | {
|
---|
180 | IndexVector first_kid, next_kid; // Linked list of children
|
---|
181 | Index postnum;
|
---|
182 | // Allocate storage for working arrays and results
|
---|
183 | first_kid.resize(n+1);
|
---|
184 | next_kid.setZero(n+1);
|
---|
185 | post.setZero(n+1);
|
---|
186 |
|
---|
187 | // Set up structure describing children
|
---|
188 | Index v, dad;
|
---|
189 | first_kid.setConstant(-1);
|
---|
190 | for (v = n-1; v >= 0; v--)
|
---|
191 | {
|
---|
192 | 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
|
---|