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 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
    137