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
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>
pruneL(const Index jcol,const IndexVector & perm_r,const Index pivrow,const Index nseg,const IndexVector & segrep,BlockIndexVector repfnz,IndexVector & xprune,GlobalLU_t & glu)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