$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 00011 /* 00012 00013 * NOTE: This file is the modified version of sp_coletree.c file in SuperLU 00014 00015 * -- SuperLU routine (version 3.1) -- 00016 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 00017 * and Lawrence Berkeley National Lab. 00018 * August 1, 2008 00019 * 00020 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 00021 * 00022 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 00023 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 00024 * 00025 * Permission is hereby granted to use or copy this program for any 00026 * purpose, provided the above notices are retained on all copies. 00027 * Permission to modify the code and to distribute modified code is 00028 * granted, provided the above notices are retained, and a notice that 00029 * the code was modified is included with the above copyright notice. 00030 */ 00031 #ifndef SPARSE_COLETREE_H 00032 #define SPARSE_COLETREE_H 00033 00034 namespace Eigen { 00035 00036 namespace internal { 00037 00039 template<typename Index, typename IndexVector> 00040 Index etree_find (Index i, IndexVector& pp) 00041 { 00042 Index p = pp(i); // Parent 00043 Index gp = pp(p); // Grand parent 00044 while (gp != p) 00045 { 00046 pp(i) = gp; // Parent pointer on find path is changed to former grand parent 00047 i = gp; 00048 p = pp(i); 00049 gp = pp(p); 00050 } 00051 return p; 00052 } 00053 00060 template <typename MatrixType, typename IndexVector> 00061 int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::Index *perm=0) 00062 { 00063 typedef typename MatrixType::Index Index; 00064 Index nc = mat.cols(); // Number of columns 00065 Index m = mat.rows(); 00066 Index diagSize = (std::min)(nc,m); 00067 IndexVector root(nc); // root of subtree of etree 00068 root.setZero(); 00069 IndexVector pp(nc); // disjoint sets 00070 pp.setZero(); // Initialize disjoint sets 00071 parent.resize(mat.cols()); 00072 //Compute first nonzero column in each row 00073 Index row,col; 00074 firstRowElt.resize(m); 00075 firstRowElt.setConstant(nc); 00076 firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1); 00077 bool found_diag; 00078 for (col = 0; col < nc; col++) 00079 { 00080 Index pcol = col; 00081 if(perm) pcol = perm[col]; 00082 for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it) 00083 { 00084 row = it.row(); 00085 firstRowElt(row) = (std::min)(firstRowElt(row), col); 00086 } 00087 } 00088 /* Compute etree by Liu's algorithm for symmetric matrices, 00089 except use (firstRowElt[r],c) in place of an edge (r,c) of A. 00090 Thus each row clique in A'*A is replaced by a star 00091 centered at its first vertex, which has the same fill. */ 00092 Index rset, cset, rroot; 00093 for (col = 0; col < nc; col++) 00094 { 00095 found_diag = col>=m; 00096 pp(col) = col; 00097 cset = col; 00098 root(cset) = col; 00099 parent(col) = nc; 00100 /* The diagonal element is treated here even if it does not exist in the matrix 00101 * hence the loop is executed once more */ 00102 Index pcol = col; 00103 if(perm) pcol = perm[col]; 00104 for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it) 00105 { // A sequence of interleaved find and union is performed 00106 Index i = col; 00107 if(it) i = it.index(); 00108 if (i == col) found_diag = true; 00109 00110 row = firstRowElt(i); 00111 if (row >= col) continue; 00112 rset = internal::etree_find(row, pp); // Find the name of the set containing row 00113 rroot = root(rset); 00114 if (rroot != col) 00115 { 00116 parent(rroot) = col; 00117 pp(cset) = rset; 00118 cset = rset; 00119 root(cset) = col; 00120 } 00121 } 00122 } 00123 return 0; 00124 } 00125 00130 template <typename Index, typename IndexVector> 00131 void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum) 00132 { 00133 Index current = n, first, next; 00134 while (postnum != n) 00135 { 00136 // No kid for the current node 00137 first = first_kid(current); 00138 00139 // no kid for the current node 00140 if (first == -1) 00141 { 00142 // Numbering this node because it has no kid 00143 post(current) = postnum++; 00144 00145 // looking for the next kid 00146 next = next_kid(current); 00147 while (next == -1) 00148 { 00149 // No more kids : back to the parent node 00150 current = parent(current); 00151 // numbering the parent node 00152 post(current) = postnum++; 00153 00154 // Get the next kid 00155 next = next_kid(current); 00156 } 00157 // stopping criterion 00158 if (postnum == n+1) return; 00159 00160 // Updating current node 00161 current = next; 00162 } 00163 else 00164 { 00165 current = first; 00166 } 00167 } 00168 } 00169 00170 00177 template <typename Index, typename IndexVector> 00178 void treePostorder(Index n, IndexVector& parent, IndexVector& post) 00179 { 00180 IndexVector first_kid, next_kid; // Linked list of children 00181 Index postnum; 00182 // Allocate storage for working arrays and results 00183 first_kid.resize(n+1); 00184 next_kid.setZero(n+1); 00185 post.setZero(n+1); 00186 00187 // Set up structure describing children 00188 Index v, dad; 00189 first_kid.setConstant(-1); 00190 for (v = n-1; v >= 0; v--) 00191 { 00192 dad = parent(v); 00193 next_kid(v) = first_kid(dad); 00194 first_kid(dad) = v; 00195 } 00196 00197 // Depth-first search from dummy root vertex #n 00198 postnum = 0; 00199 internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum); 00200 } 00201 00202 } // end namespace internal 00203 00204 } // end namespace Eigen 00205 00206 #endif // SPARSE_COLETREE_H