SparseLU_pruneL.h
Go to the documentation of this file.
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 
52 template <typename Scalar, typename StorageIndex>
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
Eigen::internal::LU_GlobalLU_t::lsub
IndexVector lsub
Definition: SparseLU_Structs.h:82
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::internal::LU_GlobalLU_t::xsup
IndexVector xsup
Definition: SparseLU_Structs.h:79
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
Eigen::internal::LU_GlobalLU_t::xlusup
IndexVector xlusup
Definition: SparseLU_Structs.h:83
Eigen::internal::SparseLUImpl::pruneL
void pruneL(const Index jcol, const IndexVector &perm_r, const Index pivrow, const Index nseg, const IndexVector &segrep, BlockIndexVector repfnz, IndexVector &xprune, GlobalLU_t &glu)
Prunes the L-structure.
Definition: SparseLU_pruneL.h:53
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
Eigen::internal::LU_GlobalLU_t::supno
IndexVector supno
Definition: SparseLU_Structs.h:80
Eigen::internal::emptyIdxLU
@ emptyIdxLU
Definition: SparseLU_Memory.h:38
Eigen::internal::LU_GlobalLU_t::xlsub
IndexVector xlsub
Definition: SparseLU_Structs.h:84
Eigen::internal::LU_GlobalLU_t::lusup
ScalarVector lusup
Definition: SparseLU_Structs.h:81
Eigen::Matrix< StorageIndex, Dynamic, 1 >
internal
Definition: BandTriangularSolver.h:13
Eigen::internal::LU_GlobalLU_t
Definition: SparseLU_Structs.h:77
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:04:33