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/* This file is a modified version of heap_relax_snode.c file in SuperLU 11 * -- SuperLU routine (version 3.0) -- 12 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 13 * and Lawrence Berkeley National Lab. 14 * October 15, 2003 15 * 16 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 17 * 18 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 19 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 20 * 21 * Permission is hereby granted to use or copy this program for any 22 * purpose, provided the above notices are retained on all copies. 23 * Permission to modify the code and to distribute modified code is 24 * granted, provided the above notices are retained, and a notice that 25 * the code was modified is included with the above copyright notice. 26 */ 27 28#ifndef SPARSELU_HEAP_RELAX_SNODE_H 29#define SPARSELU_HEAP_RELAX_SNODE_H 30 31namespace Eigen { 32namespace internal { 33 34/** 35 * \brief Identify the initial relaxed supernodes 36 * 37 * This routine applied to a symmetric elimination tree. 38 * It assumes that the matrix has been reordered according to the postorder of the etree 39 * \param n The number of columns 40 * \param et elimination tree 41 * \param relax_columns Maximum number of columns allowed in a relaxed snode 42 * \param descendants Number of descendants of each node in the etree 43 * \param relax_end last column in a supernode 44 */ 45template <typename Scalar, typename Index> 46void SparseLUImpl<Scalar,Index>::heap_relax_snode (const Index n, IndexVector& et, const Index relax_columns, IndexVector& descendants, IndexVector& relax_end) 47{ 48 49 // The etree may not be postordered, but its heap ordered 50 IndexVector post; 51 internal::treePostorder(n, et, post); // Post order etree 52 IndexVector inv_post(n+1); 53 Index i; 54 for (i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()??? 55 56 // Renumber etree in postorder 57 IndexVector iwork(n); 58 IndexVector et_save(n+1); 59 for (i = 0; i < n; ++i) 60 { 61 iwork(post(i)) = post(et(i)); 62 } 63 et_save = et; // Save the original etree 64 et = iwork; 65 66 // compute the number of descendants of each node in the etree 67 relax_end.setConstant(emptyIdxLU); 68 Index j, parent; 69 descendants.setZero(); 70 for (j = 0; j < n; j++) 71 { 72 parent = et(j); 73 if (parent != n) // not the dummy root 74 descendants(parent) += descendants(j) + 1; 75 } 76 // Identify the relaxed supernodes by postorder traversal of the etree 77 Index snode_start; // beginning of a snode 78 Index k; 79 Index nsuper_et_post = 0; // Number of relaxed snodes in postordered etree 80 Index nsuper_et = 0; // Number of relaxed snodes in the original etree 81 Index l; 82 for (j = 0; j < n; ) 83 { 84 parent = et(j); 85 snode_start = j; 86 while ( parent != n && descendants(parent) < relax_columns ) 87 { 88 j = parent; 89 parent = et(j); 90 } 91 // Found a supernode in postordered etree, j is the last column 92 ++nsuper_et_post; 93 k = n; 94 for (i = snode_start; i <= j; ++i) 95 k = (std::min)(k, inv_post(i)); 96 l = inv_post(j); 97 if ( (l - k) == (j - snode_start) ) // Same number of columns in the snode 98 { 99 // This is also a supernode in the original etree 100 relax_end(k) = l; // Record last column 101 ++nsuper_et; 102 } 103 else 104 { 105 for (i = snode_start; i <= j; ++i) 106 { 107 l = inv_post(i); 108 if (descendants(i) == 0) 109 { 110 relax_end(l) = l; 111 ++nsuper_et; 112 } 113 } 114 } 115 j++; 116 // Search for a new leaf 117 while (descendants(j) != 0 && j < n) j++; 118 } // End postorder traversal of the etree 119 120 // Recover the original etree 121 et = et_save; 122} 123 124} // end namespace internal 125 126} // end namespace Eigen 127#endif // SPARSELU_HEAP_RELAX_SNODE_H 128