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
33namespace Eigen {
34namespace 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 */
52template <typename Scalar, typename StorageIndex>
53void 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