src/SparseCore/SparseColEtree.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010
0011 /*
0012
0013  * NOTE: This file is the modified version of sp_coletree.c file in SuperLU
0014
0015  * -- SuperLU routine (version 3.1) --
0016  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0017  * and Lawrence Berkeley National Lab.
0018  * August 1, 2008
0019  *
0020  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0021  *
0022  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0023  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0024  *
0025  * Permission is hereby granted to use or copy this program for any
0026  * purpose, provided the above notices are retained on all copies.
0027  * Permission to modify the code and to distribute modified code is
0028  * granted, provided the above notices are retained, and a notice that
0029  * the code was modified is included with the above copyright notice.
0030  */
0031 #ifndef SPARSE_COLETREE_H
0032 #define SPARSE_COLETREE_H
0033
0034 namespace Eigen {
0035
0036 namespace internal {
0037
0038 /** Find the root of the tree/set containing the vertex i : Use Path halving */
0039 template<typename Index, typename IndexVector>
0040 Index etree_find (Index i, IndexVector& pp)
0041 {
0042   Index p = pp(i); // Parent
0043   Index gp = pp(p); // Grand parent
0044   while (gp != p)
0045   {
0046     pp(i) = gp; // Parent pointer on find path is changed to former grand parent
0047     i = gp;
0048     p = pp(i);
0049     gp = pp(p);
0050   }
0051   return p;
0052 }
0053
0054 /** Compute the column elimination tree of a sparse matrix
0055   * \param mat The matrix in column-major format.
0056   * \param parent The elimination tree
0057   * \param firstRowElt The column index of the first element in each row
0058   * \param perm The permutation to apply to the column of \b mat
0059   */
0060 template <typename MatrixType, typename IndexVector>
0061 int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::StorageIndex *perm=0)
0062 {
0063   typedef typename MatrixType::StorageIndex StorageIndex;
0064   StorageIndex nc = convert_index<StorageIndex>(mat.cols()); // Number of columns
0065   StorageIndex m = convert_index<StorageIndex>(mat.rows());
0066   StorageIndex diagSize = (std::min)(nc,m);
0067   IndexVector root(nc); // root of subtree of etree
0068   root.setZero();
0069   IndexVector pp(nc); // disjoint sets
0070   pp.setZero(); // Initialize disjoint sets
0071   parent.resize(mat.cols());
0072   //Compute first nonzero column in each row
0073   firstRowElt.resize(m);
0074   firstRowElt.setConstant(nc);
0075   firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1);
0076   bool found_diag;
0077   for (StorageIndex col = 0; col < nc; col++)
0078   {
0079     StorageIndex pcol = col;
0080     if(perm) pcol  = perm[col];
0081     for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it)
0082     {
0083       Index row = it.row();
0084       firstRowElt(row) = (std::min)(firstRowElt(row), col);
0085     }
0086   }
0087   /* Compute etree by Liu's algorithm for symmetric matrices,
0088           except use (firstRowElt[r],c) in place of an edge (r,c) of A.
0089     Thus each row clique in A'*A is replaced by a star
0090     centered at its first vertex, which has the same fill. */
0091   StorageIndex rset, cset, rroot;
0092   for (StorageIndex col = 0; col < nc; col++)
0093   {
0094     found_diag = col>=m;
0095     pp(col) = col;
0096     cset = col;
0097     root(cset) = col;
0098     parent(col) = nc;
0099     /* The diagonal element is treated here even if it does not exist in the matrix
0100      * hence the loop is executed once more */
0101     StorageIndex pcol = col;
0102     if(perm) pcol  = perm[col];
0103     for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it)
0104     { //  A sequence of interleaved find and union is performed
0105       Index i = col;
0106       if(it) i = it.index();
0107       if (i == col) found_diag = true;
0108
0109       StorageIndex row = firstRowElt(i);
0110       if (row >= col) continue;
0111       rset = internal::etree_find(row, pp); // Find the name of the set containing row
0112       rroot = root(rset);
0113       if (rroot != col)
0114       {
0115         parent(rroot) = col;
0116         pp(cset) = rset;
0117         cset = rset;
0118         root(cset) = col;
0119       }
0120     }
0121   }
0122   return 0;
0123 }
0124
0125 /**
0126   * Depth-first search from vertex n.  No recursion.
0127   * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
0128 */
0129 template <typename IndexVector>
0130 void nr_etdfs (typename IndexVector::Scalar n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, typename IndexVector::Scalar postnum)
0131 {
0132   typedef typename IndexVector::Scalar StorageIndex;
0133   StorageIndex current = n, first, next;
0134   while (postnum != n)
0135   {
0136     // No kid for the current node
0137     first = first_kid(current);
0138
0139     // no kid for the current node
0140     if (first == -1)
0141     {
0142       // Numbering this node because it has no kid
0143       post(current) = postnum++;
0144
0145       // looking for the next kid
0146       next = next_kid(current);
0147       while (next == -1)
0148       {
0149         // No more kids : back to the parent node
0150         current = parent(current);
0151         // numbering the parent node
0152         post(current) = postnum++;
0153
0154         // Get the next kid
0155         next = next_kid(current);
0156       }
0157       // stopping criterion
0158       if (postnum == n+1) return;
0159
0160       // Updating current node
0161       current = next;
0162     }
0163     else
0164     {
0165       current = first;
0166     }
0167   }
0168 }
0169
0170
0171 /**
0172   * \brief Post order a tree
0173   * \param n the number of nodes
0174   * \param parent Input tree
0175   * \param post postordered tree
0176   */
0177 template <typename IndexVector>
0178 void treePostorder(typename IndexVector::Scalar n, IndexVector& parent, IndexVector& post)
0179 {
0180   typedef typename IndexVector::Scalar StorageIndex;
0181   IndexVector first_kid, next_kid; // Linked list of children
0182   StorageIndex postnum;
0183   // Allocate storage for working arrays and results
0184   first_kid.resize(n+1);
0185   next_kid.setZero(n+1);
0186   post.setZero(n+1);
0187
0188   // Set up structure describing children
0189   first_kid.setConstant(-1);
0190   for (StorageIndex v = n-1; v >= 0; v--)
0191   {
0192     StorageIndex dad = parent(v);
0193     next_kid(v) = first_kid(dad);
0194     first_kid(dad) = v;
0195   }
0196
0197   // Depth-first search from dummy root vertex #n
0198   postnum = 0;
0199   internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
0200 }
0201
0202 } // end namespace internal
0203
0204 } // end namespace Eigen
0205
0206 #endif // SPARSE_COLETREE_H