1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 31 namespace Eigen { 32 namespace 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 */ 45 template <typename Scalar, typename StorageIndex> 46 void SparseLUImpl<Scalar,StorageIndex>::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(StorageIndex(n), et, post); // Post order etree 52 IndexVector inv_post(n+1); 53 for (StorageIndex i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()??? 54 55 // Renumber etree in postorder 56 IndexVector iwork(n); 57 IndexVector et_save(n+1); 58 for (Index i = 0; i < n; ++i) 59 { 60 iwork(post(i)) = post(et(i)); 61 } 62 et_save = et; // Save the original etree 63 et = iwork; 64 65 // compute the number of descendants of each node in the etree 66 relax_end.setConstant(emptyIdxLU); 67 Index j, parent; 68 descendants.setZero(); 69 for (j = 0; j < n; j++) 70 { 71 parent = et(j); 72 if (parent != n) // not the dummy root 73 descendants(parent) += descendants(j) + 1; 74 } 75 // Identify the relaxed supernodes by postorder traversal of the etree 76 Index snode_start; // beginning of a snode 77 StorageIndex k; 78 Index nsuper_et_post = 0; // Number of relaxed snodes in postordered etree 79 Index nsuper_et = 0; // Number of relaxed snodes in the original etree 80 StorageIndex l; 81 for (j = 0; j < n; ) 82 { 83 parent = et(j); 84 snode_start = j; 85 while ( parent != n && descendants(parent) < relax_columns ) 86 { 87 j = parent; 88 parent = et(j); 89 } 90 // Found a supernode in postordered etree, j is the last column 91 ++nsuper_et_post; 92 k = StorageIndex(n); 93 for (Index i = snode_start; i <= j; ++i) 94 k = (std::min)(k, inv_post(i)); 95 l = inv_post(j); 96 if ( (l - k) == (j - snode_start) ) // Same number of columns in the snode 97 { 98 // This is also a supernode in the original etree 99 relax_end(k) = l; // Record last column 100 ++nsuper_et; 101 } 102 else 103 { 104 for (Index i = snode_start; i <= j; ++i) 105 { 106 l = inv_post(i); 107 if (descendants(i) == 0) 108 { 109 relax_end(l) = l; 110 ++nsuper_et; 111 } 112 } 113 } 114 j++; 115 // Search for a new leaf 116 while (descendants(j) != 0 && j < n) j++; 117 } // End postorder traversal of the etree 118 119 // Recover the original etree 120 et = et_save; 121 } 122 123 } // end namespace internal 124 125 } // end namespace Eigen 126 #endif // SPARSELU_HEAP_RELAX_SNODE_H 127
