1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 Cdric 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 207
