src/SparseLU/SparseLU_pruneL.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 /*
0011
0012  * NOTE: This file is the modified version of [s,d,c,z]pruneL.c file in SuperLU
0013
0014  * -- SuperLU routine (version 2.0) --
0015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0016  * and Lawrence Berkeley National Lab.
0017  * November 15, 1997
0018  *
0019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0020  *
0021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0023  *
0024  * Permission is hereby granted to use or copy this program for any
0025  * purpose, provided the above notices are retained on all copies.
0026  * Permission to modify the code and to distribute modified code is
0027  * granted, provided the above notices are retained, and a notice that
0028  * the code was modified is included with the above copyright notice.
0029  */
0030 #ifndef SPARSELU_PRUNEL_H
0031 #define SPARSELU_PRUNEL_H
0032
0033 namespace RivetEigen {
0034 namespace internal {
0035
0036 /**
0037  * \brief Prunes the L-structure.
0038  *
0039  * It prunes the L-structure  of supernodes whose L-structure contains the current pivot row "pivrow"
0040  *
0041  *
0042  * \param jcol The current column of L
0043  * \param[in] perm_r Row permutation
0044  * \param[out] pivrow  The pivot row
0045  * \param nseg Number of segments
0046  * \param segrep
0047  * \param repfnz
0048  * \param[out] xprune
0049  * \param glu Global LU data
0050  *
0051  */
0052 template <typename Scalar, typename StorageIndex>
0053 void SparseLUImpl<Scalar,StorageIndex>::pruneL(const Index jcol, const IndexVector& perm_r, const Index pivrow, const Index nseg,
0054                                                const IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, GlobalLU_t& glu)
0055 {
0056   // For each supernode-rep irep in U(*,j]
0057   Index jsupno = glu.supno(jcol);
0058   Index i,irep,irep1;
0059   bool movnum, do_prune = false;
0060   Index kmin = 0, kmax = 0, minloc, maxloc,krow;
0061   for (i = 0; i < nseg; i++)
0062   {
0063     irep = segrep(i);
0064     irep1 = irep + 1;
0065     do_prune = false;
0066
0067     // Don't prune with a zero U-segment
0068     if (repfnz(irep) == emptyIdxLU) continue;
0069
0070     // If a snode overlaps with the next panel, then the U-segment
0071     // is fragmented into two parts -- irep and irep1. We should let
0072     // pruning occur at the rep-column in irep1s snode.
0073     if (glu.supno(irep) == glu.supno(irep1) ) continue; // don't prune
0074
0075     // If it has not been pruned & it has a nonz in row L(pivrow,i)
0076     if (glu.supno(irep) != jsupno )
0077     {
0078       if ( xprune (irep) >= glu.xlsub(irep1) )
0079       {
0080         kmin = glu.xlsub(irep);
0081         kmax = glu.xlsub(irep1) - 1;
0082         for (krow = kmin; krow <= kmax; krow++)
0083         {
0084           if (glu.lsub(krow) == pivrow)
0085           {
0086             do_prune = true;
0087             break;
0088           }
0089         }
0090       }
0091
0092       if (do_prune)
0093       {
0094         // do a quicksort-type partition
0095         // movnum=true means that the num values have to be exchanged
0096         movnum = false;
0097         if (irep == glu.xsup(glu.supno(irep)) ) // Snode of size 1
0098           movnum = true;
0099
0100         while (kmin <= kmax)
0101         {
0102           if (perm_r(glu.lsub(kmax)) == emptyIdxLU)
0103             kmax--;
0104           else if ( perm_r(glu.lsub(kmin)) != emptyIdxLU)
0105             kmin++;
0106           else
0107           {
0108             // kmin below pivrow (not yet pivoted), and kmax
0109             // above pivrow: interchange the two suscripts
0110             std::swap(glu.lsub(kmin), glu.lsub(kmax));
0111
0112             // If the supernode has only one column, then we
0113             // only keep one set of subscripts. For any subscript
0114             // intercnahge performed, similar interchange must be
0115             // done on the numerical values.
0116             if (movnum)
0117             {
0118               minloc = glu.xlusup(irep) + ( kmin - glu.xlsub(irep) );
0119               maxloc = glu.xlusup(irep) + ( kmax - glu.xlsub(irep) );
0120               std::swap(glu.lusup(minloc), glu.lusup(maxloc));
0121             }
0122             kmin++;
0123             kmax--;
0124           }
0125         } // end while
0126
0127         xprune(irep) = StorageIndex(kmin);  //Pruning
0128       } // end if do_prune
0129     } // end pruning
0130   } // End for each U-segment
0131 }
0132
0133 } // end namespace internal
0134 } // end namespace RivetEigen
0135
0136 #endif // SPARSELU_PRUNEL_H