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