SparseColEtree.h
Go to the documentation of this file.
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 
34 namespace Eigen {
35 
36 namespace internal {
37 
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 
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  IndexVector root(nc); // root of subtree of etree
67  root.setZero();
68  IndexVector pp(nc); // disjoint sets
69  pp.setZero(); // Initialize disjoint sets
70  parent.resize(mat.cols());
71  //Compute first nonzero column in each row
72  Index row,col;
73  firstRowElt.resize(m);
74  firstRowElt.setConstant(nc);
75  firstRowElt.segment(0, nc).setLinSpaced(nc, 0, nc-1);
76  bool found_diag;
77  for (col = 0; col < nc; col++)
78  {
79  Index pcol = col;
80  if(perm) pcol = perm[col];
81  for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
82  {
83  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  Index rset, cset, rroot;
92  for (col = 0; col < nc; col++)
93  {
94  found_diag = false;
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  Index 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  row = firstRowElt(i);
109  if (row >= col) continue;
110  rset = internal::etree_find(row, pp); // Find the name of the set containing row
111  rroot = root(rset);
112  if (rroot != col)
113  {
114  parent(rroot) = col;
115  pp(cset) = rset;
116  cset = rset;
117  root(cset) = col;
118  }
119  }
120  }
121  return 0;
122 }
123 
128 template <typename Index, typename IndexVector>
129 void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum)
130 {
131  Index current = n, first, next;
132  while (postnum != n)
133  {
134  // No kid for the current node
135  first = first_kid(current);
136 
137  // no kid for the current node
138  if (first == -1)
139  {
140  // Numbering this node because it has no kid
141  post(current) = postnum++;
142 
143  // looking for the next kid
144  next = next_kid(current);
145  while (next == -1)
146  {
147  // No more kids : back to the parent node
148  current = parent(current);
149  // numbering the parent node
150  post(current) = postnum++;
151 
152  // Get the next kid
153  next = next_kid(current);
154  }
155  // stopping criterion
156  if (postnum == n+1) return;
157 
158  // Updating current node
159  current = next;
160  }
161  else
162  {
163  current = first;
164  }
165  }
166 }
167 
168 
175 template <typename Index, typename IndexVector>
176 void treePostorder(Index n, IndexVector& parent, IndexVector& post)
177 {
178  IndexVector first_kid, next_kid; // Linked list of children
179  Index postnum;
180  // Allocate storage for working arrays and results
181  first_kid.resize(n+1);
182  next_kid.setZero(n+1);
183  post.setZero(n+1);
184 
185  // Set up structure describing children
186  Index v, dad;
187  first_kid.setConstant(-1);
188  for (v = n-1; v >= 0; v--)
189  {
190  dad = parent(v);
191  next_kid(v) = first_kid(dad);
192  first_kid(dad) = v;
193  }
194 
195  // Depth-first search from dummy root vertex #n
196  postnum = 0;
197  internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
198 }
199 
200 } // end namespace internal
201 
202 } // end namespace Eigen
203 
204 #endif // SPARSE_COLETREE_H
Index etree_find(Index i, IndexVector &pp)
void nr_etdfs(Index n, IndexVector &parent, IndexVector &first_kid, IndexVector &next_kid, IndexVector &post, Index postnum)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::Index *perm=0)
#define v
RowXpr row(Index i)
Definition: BlockMethods.h:725
Expression next(const Expression &arg)
ColXpr col(Index i)
Definition: BlockMethods.h:708
void treePostorder(Index n, IndexVector &parent, IndexVector &post)
Post order a tree.


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:07