Home | History | Annotate | Download | only in SparseLU 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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 Index>
     53 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)
     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