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 Index>
53void 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