[136] | 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
|
---|