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/* 11 12 * NOTE: This file is the modified version of [s,d,c,z]pruneL.c file in SuperLU 13 14 * -- SuperLU routine (version 2.0) -- 15 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 16 * and Lawrence Berkeley National Lab. 17 * November 15, 1997 18 * 19 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 20 * 21 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 22 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 23 * 24 * Permission is hereby granted to use or copy this program for any 25 * purpose, provided the above notices are retained on all copies. 26 * Permission to modify the code and to distribute modified code is 27 * granted, provided the above notices are retained, and a notice that 28 * the code was modified is included with the above copyright notice. 29 */ 30#ifndef SPARSELU_PRUNEL_H 31#define SPARSELU_PRUNEL_H 32 33namespace Eigen { 34namespace internal { 35 36/** 37 * \brief Prunes the L-structure. 38 * 39 * It prunes the L-structure of supernodes whose L-structure contains the current pivot row "pivrow" 40 * 41 * 42 * \param jcol The current column of L 43 * \param[in] perm_r Row permutation 44 * \param[out] pivrow The pivot row 45 * \param nseg Number of segments 46 * \param segrep 47 * \param repfnz 48 * \param[out] xprune 49 * \param glu Global LU data 50 * 51 */ 52template <typename Scalar, typename Index> 53void SparseLUImpl<Scalar,Index>::pruneL(const Index jcol, const IndexVector& perm_r, const Index pivrow, const Index nseg, const IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, GlobalLU_t& glu) 54{ 55 // For each supernode-rep irep in U(*,j] 56 Index jsupno = glu.supno(jcol); 57 Index i,irep,irep1; 58 bool movnum, do_prune = false; 59 Index kmin = 0, kmax = 0, minloc, maxloc,krow; 60 for (i = 0; i < nseg; i++) 61 { 62 irep = segrep(i); 63 irep1 = irep + 1; 64 do_prune = false; 65 66 // Don't prune with a zero U-segment 67 if (repfnz(irep) == emptyIdxLU) continue; 68 69 // If a snode overlaps with the next panel, then the U-segment 70 // is fragmented into two parts -- irep and irep1. We should let 71 // pruning occur at the rep-column in irep1s snode. 72 if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune 73 74 // If it has not been pruned & it has a nonz in row L(pivrow,i) 75 if (glu.supno(irep) != jsupno ) 76 { 77 if ( xprune (irep) >= glu.xlsub(irep1) ) 78 { 79 kmin = glu.xlsub(irep); 80 kmax = glu.xlsub(irep1) - 1; 81 for (krow = kmin; krow <= kmax; krow++) 82 { 83 if (glu.lsub(krow) == pivrow) 84 { 85 do_prune = true; 86 break; 87 } 88 } 89 } 90 91 if (do_prune) 92 { 93 // do a quicksort-type partition 94 // movnum=true means that the num values have to be exchanged 95 movnum = false; 96 if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1 97 movnum = true; 98 99 while (kmin <= kmax) 100 { 101 if (perm_r(glu.lsub(kmax)) == emptyIdxLU) 102 kmax--; 103 else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU) 104 kmin++; 105 else 106 { 107 // kmin below pivrow (not yet pivoted), and kmax 108 // above pivrow: interchange the two suscripts 109 std::swap(glu.lsub(kmin), glu.lsub(kmax)); 110 111 // If the supernode has only one column, then we 112 // only keep one set of subscripts. For any subscript 113 // intercnahge performed, similar interchange must be 114 // done on the numerical values. 115 if (movnum) 116 { 117 minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) ); 118 maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) ); 119 std::swap(glu.lusup(minloc), glu.lusup(maxloc)); 120 } 121 kmin++; 122 kmax--; 123 } 124 } // end while 125 126 xprune(irep) = kmin; //Pruning 127 } // end if do_prune 128 } // end pruning 129 } // End for each U-segment 130} 131 132} // end namespace internal 133} // end namespace Eigen 134 135#endif // SPARSELU_PRUNEL_H 136