17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library
27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra.
37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez//
47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez//
67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla
77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed
87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*
127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE: This file is the modified version of sp_coletree.c file in SuperLU
147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * -- SuperLU routine (version 3.1) --
167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Univ. of California Berkeley, Xerox Palo Alto Research Center,
177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and Lawrence Berkeley National Lab.
187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * August 1, 2008
197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Permission is hereby granted to use or copy this program for any
267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * purpose, provided the above notices are retained on all copies.
277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Permission to modify the code and to distribute modified code is
287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * granted, provided the above notices are retained, and a notice that
297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * the code was modified is included with the above copyright notice.
307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */
317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef SPARSE_COLETREE_H
327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define SPARSE_COLETREE_H
337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen {
357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal {
377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** Find the root of the tree/set containing the vertex i : Use Path halving */
397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Index, typename IndexVector>
407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezIndex etree_find (Index i, IndexVector& pp)
417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index p = pp(i); // Parent
437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index gp = pp(p); // Grand parent
447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  while (gp != p)
457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    pp(i) = gp; // Parent pointer on find path is changed to former grand parent
477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    i = gp;
487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    p = pp(i);
497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    gp = pp(p);
507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  return p;
527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}
537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** Compute the column elimination tree of a sparse matrix
557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param mat The matrix in column-major format.
567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param parent The elimination tree
577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param firstRowElt The column index of the first element in each row
587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param perm The permutation to apply to the column of \b mat
597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  */
607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename MatrixType, typename IndexVector>
617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezint coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::Index *perm=0)
627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  typedef typename MatrixType::Index Index;
647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index nc = mat.cols(); // Number of columns
657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index m = mat.rows();
667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index diagSize = (std::min)(nc,m);
677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  IndexVector root(nc); // root of subtree of etree
687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  root.setZero();
697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  IndexVector pp(nc); // disjoint sets
707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  pp.setZero(); // Initialize disjoint sets
717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  parent.resize(mat.cols());
727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  //Compute first nonzero column in each row
737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index row,col;
747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  firstRowElt.resize(m);
757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  firstRowElt.setConstant(nc);
767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  bool found_diag;
787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  for (col = 0; col < nc; col++)
797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index pcol = col;
817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    if(perm) pcol  = perm[col];
827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      row = it.row();
857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      firstRowElt(row) = (std::min)(firstRowElt(row), col);
867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  /* Compute etree by Liu's algorithm for symmetric matrices,
897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          except use (firstRowElt[r],c) in place of an edge (r,c) of A.
907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Thus each row clique in A'*A is replaced by a star
917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    centered at its first vertex, which has the same fill. */
927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index rset, cset, rroot;
937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  for (col = 0; col < nc; col++)
947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    found_diag = col>=m;
967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    pp(col) = col;
977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    cset = col;
987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    root(cset) = col;
997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    parent(col) = nc;
1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /* The diagonal element is treated here even if it does not exist in the matrix
1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * hence the loop is executed once more */
1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index pcol = col;
1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    if(perm) pcol  = perm[col];
1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    { //  A sequence of interleaved find and union is performed
1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index i = col;
1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if(it) i = it.index();
1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (i == col) found_diag = true;
1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      row = firstRowElt(i);
1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (row >= col) continue;
1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      rset = internal::etree_find(row, pp); // Find the name of the set containing row
1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      rroot = root(rset);
1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (rroot != col)
1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      {
1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        parent(rroot) = col;
1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        pp(cset) = rset;
1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        cset = rset;
1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        root(cset) = col;
1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  return 0;
1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}
1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/**
1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * Depth-first search from vertex n.  No recursion.
1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez*/
1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename Index, typename IndexVector>
1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum)
1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index current = n, first, next;
1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  while (postnum != n)
1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    // No kid for the current node
1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    first = first_kid(current);
1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    // no kid for the current node
1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    if (first == -1)
1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      // Numbering this node because it has no kid
1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      post(current) = postnum++;
1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      // looking for the next kid
1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      next = next_kid(current);
1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      while (next == -1)
1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      {
1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // No more kids : back to the parent node
1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        current = parent(current);
1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // numbering the parent node
1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        post(current) = postnum++;
1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // Get the next kid
1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        next = next_kid(current);
1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      // stopping criterion
1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (postnum == n+1) return;
1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      // Updating current node
1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      current = next;
1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    else
1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      current = first;
1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}
1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/**
1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \brief Post order a tree
1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param n the number of nodes
1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param parent Input tree
1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \param post postordered tree
1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  */
1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename Index, typename IndexVector>
1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid treePostorder(Index n, IndexVector& parent, IndexVector& post)
1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  IndexVector first_kid, next_kid; // Linked list of children
1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index postnum;
1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  // Allocate storage for working arrays and results
1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  first_kid.resize(n+1);
1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  next_kid.setZero(n+1);
1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  post.setZero(n+1);
1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  // Set up structure describing children
1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  Index v, dad;
1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  first_kid.setConstant(-1);
1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  for (v = n-1; v >= 0; v--)
1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  {
1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    dad = parent(v);
1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    next_kid(v) = first_kid(dad);
1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    first_kid(dad) = v;
1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  }
1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  // Depth-first search from dummy root vertex #n
1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  postnum = 0;
1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}
2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal
2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen
2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // SPARSE_COLETREE_H
207