SparseColEtree.h
Go to the documentation of this file.
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   IndexVector root(nc); // root of subtree of etree 
00067   root.setZero();
00068   IndexVector pp(nc); // disjoint sets 
00069   pp.setZero(); // Initialize disjoint sets 
00070   parent.resize(mat.cols());
00071   //Compute first nonzero column in each row 
00072   Index row,col; 
00073   firstRowElt.resize(m);
00074   firstRowElt.setConstant(nc);
00075   firstRowElt.segment(0, nc).setLinSpaced(nc, 0, nc-1);
00076   bool found_diag;
00077   for (col = 0; col < nc; col++)
00078   {
00079     Index pcol = col;
00080     if(perm) pcol  = perm[col];
00081     for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
00082     { 
00083       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   Index rset, cset, rroot; 
00092   for (col = 0; col < nc; col++) 
00093   {
00094     found_diag = false;
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     Index 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       row = firstRowElt(i);
00109       if (row >= col) continue; 
00110       rset = internal::etree_find(row, pp); // Find the name of the set containing row
00111       rroot = root(rset);
00112       if (rroot != col) 
00113       {
00114         parent(rroot) = col; 
00115         pp(cset) = rset; 
00116         cset = rset; 
00117         root(cset) = col; 
00118       }
00119     }
00120   }
00121   return 0;  
00122 }
00123 
00128 template <typename Index, typename IndexVector>
00129 void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum)
00130 {
00131   Index current = n, first, next;
00132   while (postnum != n) 
00133   {
00134     // No kid for the current node
00135     first = first_kid(current);
00136     
00137     // no kid for the current node
00138     if (first == -1) 
00139     {
00140       // Numbering this node because it has no kid 
00141       post(current) = postnum++;
00142       
00143       // looking for the next kid 
00144       next = next_kid(current); 
00145       while (next == -1) 
00146       {
00147         // No more kids : back to the parent node
00148         current = parent(current); 
00149         // numbering the parent node 
00150         post(current) = postnum++;
00151         
00152         // Get the next kid 
00153         next = next_kid(current); 
00154       }
00155       // stopping criterion 
00156       if (postnum == n+1) return; 
00157       
00158       // Updating current node 
00159       current = next; 
00160     }
00161     else 
00162     {
00163       current = first; 
00164     }
00165   }
00166 }
00167 
00168 
00175 template <typename Index, typename IndexVector>
00176 void treePostorder(Index n, IndexVector& parent, IndexVector& post)
00177 {
00178   IndexVector first_kid, next_kid; // Linked list of children 
00179   Index postnum; 
00180   // Allocate storage for working arrays and results 
00181   first_kid.resize(n+1); 
00182   next_kid.setZero(n+1);
00183   post.setZero(n+1);
00184   
00185   // Set up structure describing children
00186   Index v, dad; 
00187   first_kid.setConstant(-1); 
00188   for (v = n-1; v >= 0; v--) 
00189   {
00190     dad = parent(v);
00191     next_kid(v) = first_kid(dad); 
00192     first_kid(dad) = v; 
00193   }
00194   
00195   // Depth-first search from dummy root vertex #n
00196   postnum = 0; 
00197   internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
00198 }
00199 
00200 } // end namespace internal
00201 
00202 } // end namespace Eigen
00203 
00204 #endif // SPARSE_COLETREE_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:18