SparseLU_column_dfs.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]column_dfs.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_COLUMN_DFS_H
31 #define SPARSELU_COLUMN_DFS_H
32 
33 template <typename Scalar, typename Index> class SparseLUImpl;
34 namespace Eigen {
35 
36 namespace internal {
37 
38 template<typename IndexVector, typename ScalarVector>
40 {
41  typedef typename ScalarVector::Scalar Scalar;
42  typedef typename IndexVector::Scalar Index;
44  : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu), m_luImpl(luImpl)
45  {}
46  bool update_segrep(Index /*krep*/, Index /*jj*/)
47  {
48  return true;
49  }
50  void mem_expand(IndexVector& lsub, Index& nextl, Index chmark)
51  {
52  if (nextl >= m_glu.nzlmax)
53  m_luImpl.memXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions);
54  if (chmark != (m_jcol-1)) m_jsuper_ref = emptyIdxLU;
55  }
56  enum { ExpandMem = true };
57 
58  Index m_jcol;
59  Index& m_jsuper_ref;
62 };
63 
64 
92 template <typename Scalar, typename Index>
93 Index SparseLUImpl<Scalar,Index>::column_dfs(const Index m, const Index jcol, IndexVector& perm_r, Index maxsuper, Index& nseg, BlockIndexVector lsub_col, IndexVector& segrep, BlockIndexVector repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
94 {
95 
96  Index jsuper = glu.supno(jcol);
97  Index nextl = glu.xlsub(jcol);
98  VectorBlock<IndexVector> marker2(marker, 2*m, m);
99 
100 
101  column_dfs_traits<IndexVector, ScalarVector> traits(jcol, jsuper, glu, *this);
102 
103  // For each nonzero in A(*,jcol) do dfs
104  for (Index k = 0; ((k < m) ? lsub_col[k] != emptyIdxLU : false) ; k++)
105  {
106  Index krow = lsub_col(k);
107  lsub_col(k) = emptyIdxLU;
108  Index kmark = marker2(krow);
109 
110  // krow was visited before, go to the next nonz;
111  if (kmark == jcol) continue;
112 
113  dfs_kernel(jcol, perm_r, nseg, glu.lsub, segrep, repfnz, xprune, marker2, parent,
114  xplore, glu, nextl, krow, traits);
115  } // for each nonzero ...
116 
117  Index fsupc, jptr, jm1ptr, ito, ifrom, istop;
118  Index nsuper = glu.supno(jcol);
119  Index jcolp1 = jcol + 1;
120  Index jcolm1 = jcol - 1;
121 
122  // check to see if j belongs in the same supernode as j-1
123  if ( jcol == 0 )
124  { // Do nothing for column 0
125  nsuper = glu.supno(0) = 0 ;
126  }
127  else
128  {
129  fsupc = glu.xsup(nsuper);
130  jptr = glu.xlsub(jcol); // Not yet compressed
131  jm1ptr = glu.xlsub(jcolm1);
132 
133  // Use supernodes of type T2 : see SuperLU paper
134  if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = emptyIdxLU;
135 
136  // Make sure the number of columns in a supernode doesn't
137  // exceed threshold
138  if ( (jcol - fsupc) >= maxsuper) jsuper = emptyIdxLU;
139 
140  /* If jcol starts a new supernode, reclaim storage space in
141  * glu.lsub from previous supernode. Note we only store
142  * the subscript set of the first and last columns of
143  * a supernode. (first for num values, last for pruning)
144  */
145  if (jsuper == emptyIdxLU)
146  { // starts a new supernode
147  if ( (fsupc < jcolm1-1) )
148  { // >= 3 columns in nsuper
149  ito = glu.xlsub(fsupc+1);
150  glu.xlsub(jcolm1) = ito;
151  istop = ito + jptr - jm1ptr;
152  xprune(jcolm1) = istop; // intialize xprune(jcol-1)
153  glu.xlsub(jcol) = istop;
154 
155  for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
156  glu.lsub(ito) = glu.lsub(ifrom);
157  nextl = ito; // = istop + length(jcol)
158  }
159  nsuper++;
160  glu.supno(jcol) = nsuper;
161  } // if a new supernode
162  } // end else: jcol > 0
163 
164  // Tidy up the pointers before exit
165  glu.xsup(nsuper+1) = jcolp1;
166  glu.supno(jcolp1) = nsuper;
167  xprune(jcol) = nextl; // Intialize upper bound for pruning
168  glu.xlsub(jcolp1) = nextl;
169 
170  return 0;
171 }
172 
173 } // end namespace internal
174 
175 } // end namespace Eigen
176 
177 #endif
column_dfs_traits(Index jcol, Index &jsuper, typename SparseLUImpl< Scalar, Index >::GlobalLU_t &glu, SparseLUImpl< Scalar, Index > &luImpl)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Index column_dfs(const Index m, const Index jcol, IndexVector &perm_r, Index maxsuper, Index &nseg, BlockIndexVector lsub_col, IndexVector &segrep, BlockIndexVector repfnz, IndexVector &xprune, IndexVector &marker, IndexVector &parent, IndexVector &xplore, GlobalLU_t &glu)
Performs a symbolic factorization on column jcol and decide the supernode boundary.
Expression of a fixed-size or dynamic-size sub-vector.
A matrix or vector expression mapping an existing expressions.
Definition: Ref.h:17
SparseLUImpl< Scalar, Index > & m_luImpl
void mem_expand(IndexVector &lsub, Index &nextl, Index chmark)
SparseLUImpl< Scalar, Index >::GlobalLU_t & m_glu


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:08