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