Home | History | Annotate | Download | only in SparseCore
      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