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