/usr/include/eigen3/Eigen/src/SparseCore/SparseColEtree.h is in libeigen3-dev 3.2.2-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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
|