SparseColEtree.h 6.03 KB
Newer Older
Don Gagne's avatar
Don Gagne committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.


/* 
 
 * NOTE: This file is the modified version of sp_coletree.c file in SuperLU 
 
 * -- SuperLU routine (version 3.1) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * August 1, 2008
 *
 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 *
 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 *
 * Permission is hereby granted to use or copy this program for any
 * purpose, provided the above notices are retained on all copies.
 * Permission to modify the code and to distribute modified code is
 * granted, provided the above notices are retained, and a notice that
 * the code was modified is included with the above copyright notice.
 */
#ifndef SPARSE_COLETREE_H
#define SPARSE_COLETREE_H

namespace Eigen {

namespace internal {

/** Find the root of the tree/set containing the vertex i : Use Path halving */ 
template<typename Index, typename IndexVector>
Index etree_find (Index i, IndexVector& pp)
{
  Index p = pp(i); // Parent 
  Index gp = pp(p); // Grand parent 
  while (gp != p) 
  {
    pp(i) = gp; // Parent pointer on find path is changed to former grand parent
    i = gp; 
    p = pp(i);
    gp = pp(p);
  }
  return p; 
}

/** Compute the column elimination tree of a sparse matrix
  * \param mat The matrix in column-major format. 
  * \param parent The elimination tree
  * \param firstRowElt The column index of the first element in each row
  * \param perm The permutation to apply to the column of \b mat
  */
template <typename MatrixType, typename IndexVector>
int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::Index *perm=0)
{
  typedef typename MatrixType::Index Index;
  Index nc = mat.cols(); // Number of columns 
  Index m = mat.rows();
66
  Index diagSize = (std::min)(nc,m);
Don Gagne's avatar
Don Gagne committed
67 68 69 70 71 72 73 74 75
  IndexVector root(nc); // root of subtree of etree 
  root.setZero();
  IndexVector pp(nc); // disjoint sets 
  pp.setZero(); // Initialize disjoint sets 
  parent.resize(mat.cols());
  //Compute first nonzero column in each row 
  Index row,col; 
  firstRowElt.resize(m);
  firstRowElt.setConstant(nc);
76
  firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
Don Gagne's avatar
Don Gagne committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
  bool found_diag;
  for (col = 0; col < nc; col++)
  {
    Index pcol = col;
    if(perm) pcol  = perm[col];
    for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
    { 
      row = it.row();
      firstRowElt(row) = (std::min)(firstRowElt(row), col);
    }
  }
  /* Compute etree by Liu's algorithm for symmetric matrices,
          except use (firstRowElt[r],c) in place of an edge (r,c) of A.
    Thus each row clique in A'*A is replaced by a star
    centered at its first vertex, which has the same fill. */
  Index rset, cset, rroot; 
  for (col = 0; col < nc; col++) 
  {
95
    found_diag = col>=m;
Don Gagne's avatar
Don Gagne committed
96 97 98 99 100 101 102 103 104 105 106 107 108
    pp(col) = col; 
    cset = col; 
    root(cset) = col; 
    parent(col) = nc; 
    /* The diagonal element is treated here even if it does not exist in the matrix
     * hence the loop is executed once more */ 
    Index pcol = col;
    if(perm) pcol  = perm[col];
    for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
    { //  A sequence of interleaved find and union is performed 
      Index i = col;
      if(it) i = it.index();
      if (i == col) found_diag = true;
109
      
Don Gagne's avatar
Don Gagne committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
      row = firstRowElt(i);
      if (row >= col) continue; 
      rset = internal::etree_find(row, pp); // Find the name of the set containing row
      rroot = root(rset);
      if (rroot != col) 
      {
        parent(rroot) = col; 
        pp(cset) = rset; 
        cset = rset; 
        root(cset) = col; 
      }
    }
  }
  return 0;  
}

/** 
  * Depth-first search from vertex n.  No recursion.
  * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
*/
template <typename Index, typename IndexVector>
void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum)
{
  Index current = n, first, next;
  while (postnum != n) 
  {
    // No kid for the current node
    first = first_kid(current);
    
    // no kid for the current node
    if (first == -1) 
    {
      // Numbering this node because it has no kid 
      post(current) = postnum++;
      
      // looking for the next kid 
      next = next_kid(current); 
      while (next == -1) 
      {
        // No more kids : back to the parent node
        current = parent(current); 
        // numbering the parent node 
        post(current) = postnum++;
        
        // Get the next kid 
        next = next_kid(current); 
      }
      // stopping criterion 
      if (postnum == n+1) return; 
      
      // Updating current node 
      current = next; 
    }
    else 
    {
      current = first; 
    }
  }
}


/**
  * \brief Post order a tree 
  * \param n the number of nodes
  * \param parent Input tree
  * \param post postordered tree
  */
template <typename Index, typename IndexVector>
void treePostorder(Index n, IndexVector& parent, IndexVector& post)
{
  IndexVector first_kid, next_kid; // Linked list of children 
  Index postnum; 
  // Allocate storage for working arrays and results 
  first_kid.resize(n+1); 
  next_kid.setZero(n+1);
  post.setZero(n+1);
  
  // Set up structure describing children
  Index v, dad; 
  first_kid.setConstant(-1); 
  for (v = n-1; v >= 0; v--) 
  {
    dad = parent(v);
    next_kid(v) = first_kid(dad); 
    first_kid(dad) = v; 
  }
  
  // Depth-first search from dummy root vertex #n
  postnum = 0; 
  internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
}

} // end namespace internal

} // end namespace Eigen

#endif // SPARSE_COLETREE_H