SparseLU_pruneL.h
Go to the documentation of this file.
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 Index>
00053 void 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)
00054 {
00055   // For each supernode-rep irep in U(*,j]
00056   Index jsupno = glu.supno(jcol); 
00057   Index i,irep,irep1; 
00058   bool movnum, do_prune = false; 
00059   Index kmin = 0, kmax = 0, minloc, maxloc,krow; 
00060   for (i = 0; i < nseg; i++)
00061   {
00062     irep = segrep(i); 
00063     irep1 = irep + 1; 
00064     do_prune = false; 
00065     
00066     // Don't prune with a zero U-segment 
00067     if (repfnz(irep) == emptyIdxLU) continue; 
00068     
00069     // If a snode overlaps with the next panel, then the U-segment
00070     // is fragmented into two parts -- irep and irep1. We should let 
00071     // pruning occur at the rep-column in irep1s snode. 
00072     if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune 
00073     
00074     // If it has not been pruned & it has a nonz in row L(pivrow,i)
00075     if (glu.supno(irep) != jsupno )
00076     {
00077       if ( xprune (irep) >= glu.xlsub(irep1) )
00078       {
00079         kmin = glu.xlsub(irep);
00080         kmax = glu.xlsub(irep1) - 1; 
00081         for (krow = kmin; krow <= kmax; krow++)
00082         {
00083           if (glu.lsub(krow) == pivrow) 
00084           {
00085             do_prune = true; 
00086             break; 
00087           }
00088         }
00089       }
00090       
00091       if (do_prune) 
00092       {
00093         // do a quicksort-type partition
00094         // movnum=true means that the num values have to be exchanged
00095         movnum = false; 
00096         if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1 
00097           movnum = true; 
00098         
00099         while (kmin <= kmax)
00100         {
00101           if (perm_r(glu.lsub(kmax)) == emptyIdxLU)
00102             kmax--; 
00103           else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU)
00104             kmin++;
00105           else 
00106           {
00107             // kmin below pivrow (not yet pivoted), and kmax
00108             // above pivrow: interchange the two suscripts
00109             std::swap(glu.lsub(kmin), glu.lsub(kmax)); 
00110             
00111             // If the supernode has only one column, then we 
00112             // only keep one set of subscripts. For any subscript
00113             // intercnahge performed, similar interchange must be 
00114             // done on the numerical values. 
00115             if (movnum) 
00116             {
00117               minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) ); 
00118               maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) ); 
00119               std::swap(glu.lusup(minloc), glu.lusup(maxloc)); 
00120             }
00121             kmin++;
00122             kmax--;
00123           }
00124         } // end while 
00125         
00126         xprune(irep) = kmin;  //Pruning 
00127       } // end if do_prune 
00128     } // end pruning 
00129   } // End for each U-segment
00130 }
00131 
00132 } // end namespace internal
00133 } // end namespace Eigen
00134 
00135 #endif // SPARSELU_PRUNEL_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:23