![]() |
Eigen
3.3.3
|
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::StorageIndex *perm=0) 00062 { 00063 typedef typename MatrixType::StorageIndex StorageIndex; 00064 StorageIndex nc = convert_index<StorageIndex>(mat.cols()); // Number of columns 00065 StorageIndex m = convert_index<StorageIndex>(mat.rows()); 00066 StorageIndex 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 firstRowElt.resize(m); 00074 firstRowElt.setConstant(nc); 00075 firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1); 00076 bool found_diag; 00077 for (StorageIndex col = 0; col < nc; col++) 00078 { 00079 StorageIndex pcol = col; 00080 if(perm) pcol = perm[col]; 00081 for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it) 00082 { 00083 Index row = it.row(); 00084 firstRowElt(row) = (std::min)(firstRowElt(row), col); 00085 } 00086 } 00087 /* Compute etree by Liu's algorithm for symmetric matrices, 00088 except use (firstRowElt[r],c) in place of an edge (r,c) of A. 00089 Thus each row clique in A'*A is replaced by a star 00090 centered at its first vertex, which has the same fill. */ 00091 StorageIndex rset, cset, rroot; 00092 for (StorageIndex col = 0; col < nc; col++) 00093 { 00094 found_diag = col>=m; 00095 pp(col) = col; 00096 cset = col; 00097 root(cset) = col; 00098 parent(col) = nc; 00099 /* The diagonal element is treated here even if it does not exist in the matrix 00100 * hence the loop is executed once more */ 00101 StorageIndex pcol = col; 00102 if(perm) pcol = perm[col]; 00103 for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it) 00104 { // A sequence of interleaved find and union is performed 00105 Index i = col; 00106 if(it) i = it.index(); 00107 if (i == col) found_diag = true; 00108 00109 StorageIndex row = firstRowElt(i); 00110 if (row >= col) continue; 00111 rset = internal::etree_find(row, pp); // Find the name of the set containing row 00112 rroot = root(rset); 00113 if (rroot != col) 00114 { 00115 parent(rroot) = col; 00116 pp(cset) = rset; 00117 cset = rset; 00118 root(cset) = col; 00119 } 00120 } 00121 } 00122 return 0; 00123 } 00124 00129 template <typename IndexVector> 00130 void nr_etdfs (typename IndexVector::Scalar n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, typename IndexVector::Scalar postnum) 00131 { 00132 typedef typename IndexVector::Scalar StorageIndex; 00133 StorageIndex 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 IndexVector> 00178 void treePostorder(typename IndexVector::Scalar n, IndexVector& parent, IndexVector& post) 00179 { 00180 typedef typename IndexVector::Scalar StorageIndex; 00181 IndexVector first_kid, next_kid; // Linked list of children 00182 StorageIndex postnum; 00183 // Allocate storage for working arrays and results 00184 first_kid.resize(n+1); 00185 next_kid.setZero(n+1); 00186 post.setZero(n+1); 00187 00188 // Set up structure describing children 00189 first_kid.setConstant(-1); 00190 for (StorageIndex v = n-1; v >= 0; v--) 00191 { 00192 StorageIndex 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