src/SparseLU/SparseLU_heap_relax_snode.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 /* This file is a modified version of heap_relax_snode.c file in SuperLU
0011  * -- SuperLU routine (version 3.0) --
0012  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0013  * and Lawrence Berkeley National Lab.
0014  * October 15, 2003
0015  *
0016  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0017  *
0018  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0019  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0020  *
0021  * Permission is hereby granted to use or copy this program for any
0022  * purpose, provided the above notices are retained on all copies.
0023  * Permission to modify the code and to distribute modified code is
0024  * granted, provided the above notices are retained, and a notice that
0025  * the code was modified is included with the above copyright notice.
0026  */
0027
0028 #ifndef SPARSELU_HEAP_RELAX_SNODE_H
0029 #define SPARSELU_HEAP_RELAX_SNODE_H
0030
0031 namespace Eigen {
0032 namespace internal {
0033
0034 /**
0035  * \brief Identify the initial relaxed supernodes
0036  *
0037  * This routine applied to a symmetric elimination tree.
0038  * It assumes that the matrix has been reordered according to the postorder of the etree
0039  * \param n The number of columns
0040  * \param et elimination tree
0041  * \param relax_columns Maximum number of columns allowed in a relaxed snode
0042  * \param descendants Number of descendants of each node in the etree
0043  * \param relax_end last column in a supernode
0044  */
0045 template <typename Scalar, typename StorageIndex>
0046 void SparseLUImpl<Scalar,StorageIndex>::heap_relax_snode (const Index n, IndexVector& et, const Index relax_columns, IndexVector& descendants, IndexVector& relax_end)
0047 {
0048
0049   // The etree may not be postordered, but its heap ordered
0050   IndexVector post;
0051   internal::treePostorder(StorageIndex(n), et, post); // Post order etree
0052   IndexVector inv_post(n+1);
0053   for (StorageIndex i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()???
0054
0055   // Renumber etree in postorder
0056   IndexVector iwork(n);
0057   IndexVector et_save(n+1);
0058   for (Index i = 0; i < n; ++i)
0059   {
0060     iwork(post(i)) = post(et(i));
0061   }
0062   et_save = et; // Save the original etree
0063   et = iwork;
0064
0065   // compute the number of descendants of each node in the etree
0066   relax_end.setConstant(emptyIdxLU);
0067   Index j, parent;
0068   descendants.setZero();
0069   for (j = 0; j < n; j++)
0070   {
0071     parent = et(j);
0072     if (parent != n) // not the dummy root
0073       descendants(parent) += descendants(j) + 1;
0074   }
0075   // Identify the relaxed supernodes by postorder traversal of the etree
0076   Index snode_start; // beginning of a snode
0077   StorageIndex k;
0078   Index nsuper_et_post = 0; // Number of relaxed snodes in postordered etree
0079   Index nsuper_et = 0; // Number of relaxed snodes in the original etree
0080   StorageIndex l;
0081   for (j = 0; j < n; )
0082   {
0083     parent = et(j);
0084     snode_start = j;
0085     while ( parent != n && descendants(parent) < relax_columns )
0086     {
0087       j = parent;
0088       parent = et(j);
0089     }
0090     // Found a supernode in postordered etree, j is the last column
0091     ++nsuper_et_post;
0092     k = StorageIndex(n);
0093     for (Index i = snode_start; i <= j; ++i)
0094       k = (std::min)(k, inv_post(i));
0095     l = inv_post(j);
0096     if ( (l - k) == (j - snode_start) )  // Same number of columns in the snode
0097     {
0098       // This is also a supernode in the original etree
0099       relax_end(k) = l; // Record last column
0100       ++nsuper_et;
0101     }
0102     else
0103     {
0104       for (Index i = snode_start; i <= j; ++i)
0105       {
0106         l = inv_post(i);
0107         if (descendants(i) == 0)
0108         {
0109           relax_end(l) = l;
0110           ++nsuper_et;
0111         }
0112       }
0113     }
0114     j++;
0115     // Search for a new leaf
0116     while (descendants(j) != 0 && j < n) j++;
0117   } // End postorder traversal of the etree
0118
0119   // Recover the original etree
0120   et = et_save;
0121 }
0122
0123 } // end namespace internal
0124
0125 } // end namespace Eigen
0126 #endif // SPARSELU_HEAP_RELAX_SNODE_H