Eigen  3.3.3
SparseLU_pruneL.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends