1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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
34namespace Eigen {
35
36namespace internal {
37
38/** Find the root of the tree/set containing the vertex i : Use Path halving */
39template<typename Index, typename IndexVector>
40Index 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  */
60template <typename MatrixType, typename IndexVector>
61int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::StorageIndex *perm=0)
62{
63  typedef typename MatrixType::StorageIndex StorageIndex;
64  StorageIndex nc = convert_index<StorageIndex>(mat.cols()); // Number of columns
65  StorageIndex m = convert_index<StorageIndex>(mat.rows());
66  StorageIndex 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  firstRowElt.resize(m);
74  firstRowElt.setConstant(nc);
75  firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
76  bool found_diag;
77  for (StorageIndex col = 0; col < nc; col++)
78  {
79    StorageIndex pcol = col;
80    if(perm) pcol  = perm[col];
81    for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
82    {
83      Index row = it.row();
84      firstRowElt(row) = (std::min)(firstRowElt(row), col);
85    }
86  }
87  /* Compute etree by Liu's algorithm for symmetric matrices,
88          except use (firstRowElt[r],c) in place of an edge (r,c) of A.
89    Thus each row clique in A'*A is replaced by a star
90    centered at its first vertex, which has the same fill. */
91  StorageIndex rset, cset, rroot;
92  for (StorageIndex col = 0; col < nc; col++)
93  {
94    found_diag = col>=m;
95    pp(col) = col;
96    cset = col;
97    root(cset) = col;
98    parent(col) = nc;
99    /* The diagonal element is treated here even if it does not exist in the matrix
100     * hence the loop is executed once more */
101    StorageIndex pcol = col;
102    if(perm) pcol  = perm[col];
103    for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
104    { //  A sequence of interleaved find and union is performed
105      Index i = col;
106      if(it) i = it.index();
107      if (i == col) found_diag = true;
108
109      StorageIndex row = firstRowElt(i);
110      if (row >= col) continue;
111      rset = internal::etree_find(row, pp); // Find the name of the set containing row
112      rroot = root(rset);
113      if (rroot != col)
114      {
115        parent(rroot) = col;
116        pp(cset) = rset;
117        cset = rset;
118        root(cset) = col;
119      }
120    }
121  }
122  return 0;
123}
124
125/**
126  * Depth-first search from vertex n.  No recursion.
127  * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
128*/
129template <typename IndexVector>
130void nr_etdfs (typename IndexVector::Scalar n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, typename IndexVector::Scalar postnum)
131{
132  typedef typename IndexVector::Scalar StorageIndex;
133  StorageIndex 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  */
177template <typename IndexVector>
178void treePostorder(typename IndexVector::Scalar n, IndexVector& parent, IndexVector& post)
179{
180  typedef typename IndexVector::Scalar StorageIndex;
181  IndexVector first_kid, next_kid; // Linked list of children
182  StorageIndex postnum;
183  // Allocate storage for working arrays and results
184  first_kid.resize(n+1);
185  next_kid.setZero(n+1);
186  post.setZero(n+1);
187
188  // Set up structure describing children
189  first_kid.setConstant(-1);
190  for (StorageIndex v = n-1; v >= 0; v--)
191  {
192    StorageIndex 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