![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 /* 00011 00012 * NOTE: This file is the modified version of [s,d,c,z]pruneL.c file in SuperLU 00013 00014 * -- SuperLU routine (version 2.0) -- 00015 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 00016 * and Lawrence Berkeley National Lab. 00017 * November 15, 1997 00018 * 00019 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 00020 * 00021 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 00022 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 00023 * 00024 * Permission is hereby granted to use or copy this program for any 00025 * purpose, provided the above notices are retained on all copies. 00026 * Permission to modify the code and to distribute modified code is 00027 * granted, provided the above notices are retained, and a notice that 00028 * the code was modified is included with the above copyright notice. 00029 */ 00030 #ifndef SPARSELU_PRUNEL_H 00031 #define SPARSELU_PRUNEL_H 00032 00033 namespace Eigen { 00034 namespace internal { 00035 00052 template <typename Scalar, typename StorageIndex> 00053 void SparseLUImpl<Scalar,StorageIndex>::pruneL(const Index jcol, const IndexVector& perm_r, const Index pivrow, const Index nseg, 00054 const IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, GlobalLU_t& glu) 00055 { 00056 // For each supernode-rep irep in U(*,j] 00057 Index jsupno = glu.supno(jcol); 00058 Index i,irep,irep1; 00059 bool movnum, do_prune = false; 00060 Index kmin = 0, kmax = 0, minloc, maxloc,krow; 00061 for (i = 0; i < nseg; i++) 00062 { 00063 irep = segrep(i); 00064 irep1 = irep + 1; 00065 do_prune = false; 00066 00067 // Don't prune with a zero U-segment 00068 if (repfnz(irep) == emptyIdxLU) continue; 00069 00070 // If a snode overlaps with the next panel, then the U-segment 00071 // is fragmented into two parts -- irep and irep1. We should let 00072 // pruning occur at the rep-column in irep1s snode. 00073 if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune 00074 00075 // If it has not been pruned & it has a nonz in row L(pivrow,i) 00076 if (glu.supno(irep) != jsupno ) 00077 { 00078 if ( xprune (irep) >= glu.xlsub(irep1) ) 00079 { 00080 kmin = glu.xlsub(irep); 00081 kmax = glu.xlsub(irep1) - 1; 00082 for (krow = kmin; krow <= kmax; krow++) 00083 { 00084 if (glu.lsub(krow) == pivrow) 00085 { 00086 do_prune = true; 00087 break; 00088 } 00089 } 00090 } 00091 00092 if (do_prune) 00093 { 00094 // do a quicksort-type partition 00095 // movnum=true means that the num values have to be exchanged 00096 movnum = false; 00097 if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1 00098 movnum = true; 00099 00100 while (kmin <= kmax) 00101 { 00102 if (perm_r(glu.lsub(kmax)) == emptyIdxLU) 00103 kmax--; 00104 else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU) 00105 kmin++; 00106 else 00107 { 00108 // kmin below pivrow (not yet pivoted), and kmax 00109 // above pivrow: interchange the two suscripts 00110 std::swap(glu.lsub(kmin), glu.lsub(kmax)); 00111 00112 // If the supernode has only one column, then we 00113 // only keep one set of subscripts. For any subscript 00114 // intercnahge performed, similar interchange must be 00115 // done on the numerical values. 00116 if (movnum) 00117 { 00118 minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) ); 00119 maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) ); 00120 std::swap(glu.lusup(minloc), glu.lusup(maxloc)); 00121 } 00122 kmin++; 00123 kmax--; 00124 } 00125 } // end while 00126 00127 xprune(irep) = StorageIndex(kmin); //Pruning 00128 } // end if do_prune 00129 } // end pruning 00130 } // End for each U-segment 00131 } 00132 00133 } // end namespace internal 00134 } // end namespace Eigen 00135 00136 #endif // SPARSELU_PRUNEL_H