cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

SparseLU_pruneL.h (4545B)

      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 
     33 namespace Eigen {
     34 namespace 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  */
     52 template <typename Scalar, typename StorageIndex>
     53 void SparseLUImpl<Scalar,StorageIndex>::pruneL(const Index jcol, const IndexVector& perm_r, const Index pivrow, const Index nseg,
     54                                                const IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, GlobalLU_t& glu)
     55 {
     56   // For each supernode-rep irep in U(*,j]
     57   Index jsupno = glu.supno(jcol); 
     58   Index i,irep,irep1; 
     59   bool movnum, do_prune = false; 
     60   Index kmin = 0, kmax = 0, minloc, maxloc,krow; 
     61   for (i = 0; i < nseg; i++)
     62   {
     63     irep = segrep(i); 
     64     irep1 = irep + 1; 
     65     do_prune = false; 
     66     
     67     // Don't prune with a zero U-segment 
     68     if (repfnz(irep) == emptyIdxLU) continue; 
     69     
     70     // If a snode overlaps with the next panel, then the U-segment
     71     // is fragmented into two parts -- irep and irep1. We should let 
     72     // pruning occur at the rep-column in irep1s snode. 
     73     if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune 
     74     
     75     // If it has not been pruned & it has a nonz in row L(pivrow,i)
     76     if (glu.supno(irep) != jsupno )
     77     {
     78       if ( xprune (irep) >= glu.xlsub(irep1) )
     79       {
     80         kmin = glu.xlsub(irep);
     81         kmax = glu.xlsub(irep1) - 1; 
     82         for (krow = kmin; krow <= kmax; krow++)
     83         {
     84           if (glu.lsub(krow) == pivrow) 
     85           {
     86             do_prune = true; 
     87             break; 
     88           }
     89         }
     90       }
     91       
     92       if (do_prune) 
     93       {
     94         // do a quicksort-type partition
     95         // movnum=true means that the num values have to be exchanged
     96         movnum = false; 
     97         if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1 
     98           movnum = true; 
     99         
    100         while (kmin <= kmax)
    101         {
    102           if (perm_r(glu.lsub(kmax)) == emptyIdxLU)
    103             kmax--; 
    104           else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU)
    105             kmin++;
    106           else 
    107           {
    108             // kmin below pivrow (not yet pivoted), and kmax
    109             // above pivrow: interchange the two suscripts
    110             std::swap(glu.lsub(kmin), glu.lsub(kmax)); 
    111             
    112             // If the supernode has only one column, then we 
    113             // only keep one set of subscripts. For any subscript
    114             // intercnahge performed, similar interchange must be 
    115             // done on the numerical values. 
    116             if (movnum) 
    117             {
    118               minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) ); 
    119               maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) ); 
    120               std::swap(glu.lusup(minloc), glu.lusup(maxloc)); 
    121             }
    122             kmin++;
    123             kmax--;
    124           }
    125         } // end while 
    126         
    127         xprune(irep) = StorageIndex(kmin);  //Pruning 
    128       } // end if do_prune 
    129     } // end pruning 
    130   } // End for each U-segment
    131 }
    132 
    133 } // end namespace internal
    134 } // end namespace Eigen
    135 
    136 #endif // SPARSELU_PRUNEL_H