SparseLU_Memory.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]memory.c files in SuperLU
13 
14  * -- SuperLU routine (version 3.1) --
15  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
16  * and Lawrence Berkeley National Lab.
17  * August 1, 2008
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 
31 #ifndef EIGEN_SPARSELU_MEMORY
32 #define EIGEN_SPARSELU_MEMORY
33 
34 namespace Eigen {
35 namespace internal {
36 
37 enum { LUNoMarker = 3 };
38 enum {emptyIdxLU = -1};
40 {
41  return (std::max)(m, (t+b)*w);
42 }
43 
44 template< typename Scalar>
46 {
47  return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar);
48 }
49 
50 
51 
52 
61 template <typename Scalar, typename StorageIndex>
62 template <typename VectorType>
63 Index SparseLUImpl<Scalar,StorageIndex>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions)
64 {
65 
66  float alpha = 1.5; // Ratio of the memory increase
67  Index new_len; // New size of the allocated memory
68 
69  if(num_expansions == 0 || keep_prev)
70  new_len = length ; // First time allocate requested
71  else
72  new_len = (std::max)(length+1,Index(alpha * length));
73 
74  VectorType old_vec; // Temporary vector to hold the previous values
75  if (nbElts > 0 )
76  old_vec = vec.segment(0,nbElts);
77 
78  //Allocate or expand the current vector
79 #ifdef EIGEN_EXCEPTIONS
80  try
81 #endif
82  {
83  vec.resize(new_len);
84  }
85 #ifdef EIGEN_EXCEPTIONS
86  catch(std::bad_alloc& )
87 #else
88  if(!vec.size())
89 #endif
90  {
91  if (!num_expansions)
92  {
93  // First time to allocate from LUMemInit()
94  // Let LUMemInit() deals with it.
95  return -1;
96  }
97  if (keep_prev)
98  {
99  // In this case, the memory length should not not be reduced
100  return new_len;
101  }
102  else
103  {
104  // Reduce the size and increase again
105  Index tries = 0; // Number of attempts
106  do
107  {
108  alpha = (alpha + 1)/2;
109  new_len = (std::max)(length+1,Index(alpha * length));
110 #ifdef EIGEN_EXCEPTIONS
111  try
112 #endif
113  {
114  vec.resize(new_len);
115  }
116 #ifdef EIGEN_EXCEPTIONS
117  catch(std::bad_alloc& )
118 #else
119  if (!vec.size())
120 #endif
121  {
122  tries += 1;
123  if ( tries > 10) return new_len;
124  }
125  } while (!vec.size());
126  }
127  }
128  //Copy the previous values to the newly allocated space
129  if (nbElts > 0)
130  vec.segment(0, nbElts) = old_vec;
131 
132 
133  length = new_len;
134  if(num_expansions) ++num_expansions;
135  return 0;
136 }
137 
150 template <typename Scalar, typename StorageIndex>
152 {
153  Index& num_expansions = glu.num_expansions; //No memory expansions so far
154  num_expansions = 0;
155  glu.nzumax = glu.nzlumax = (std::min)(fillratio * (annz+1) / n, m) * n; // estimated number of nonzeros in U
156  glu.nzlmax = (std::max)(Index(4), fillratio) * (annz+1) / 4; // estimated nnz in L factor
157  // Return the estimated size to the user if necessary
158  Index tempSpace;
159  tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar);
160  if (lwork == emptyIdxLU)
161  {
162  Index estimated_size;
163  estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace
164  + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) * sizeof(Scalar) + n;
165  return estimated_size;
166  }
167 
168  // Setup the required space
169 
170  // First allocate Integer pointers for L\U factors
171  glu.xsup.resize(n+1);
172  glu.supno.resize(n+1);
173  glu.xlsub.resize(n+1);
174  glu.xlusup.resize(n+1);
175  glu.xusub.resize(n+1);
176 
177  // Reserve memory for L/U factors
178  do
179  {
180  if( (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0)
181  || (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0)
182  || (expand<IndexVector> (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0)
183  || (expand<IndexVector> (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) )
184  {
185  //Reduce the estimated size and retry
186  glu.nzlumax /= 2;
187  glu.nzumax /= 2;
188  glu.nzlmax /= 2;
189  if (glu.nzlumax < annz ) return glu.nzlumax;
190  }
191  } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
192 
193  ++num_expansions;
194  return 0;
195 
196 } // end LuMemInit
197 
207 template <typename Scalar, typename StorageIndex>
208 template <typename VectorType>
209 Index SparseLUImpl<Scalar,StorageIndex>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions)
210 {
211  Index failed_size;
212  if (memtype == USUB)
213  failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions);
214  else
215  failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions);
216 
217  if (failed_size)
218  return failed_size;
219 
220  return 0 ;
221 }
222 
223 } // end namespace internal
224 
225 } // end namespace Eigen
226 #endif // EIGEN_SPARSELU_MEMORY
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
Index LUTempSpace(Index &m, Index &w)
Scalar * b
Definition: benchVecAdd.cpp:17
#define min(a, b)
Definition: datatypes.h:19
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
RealScalar alpha
Index LUnumTempV(Index &m, Index &w, Index &t, Index &b)
Index memXpand(VectorType &vec, Index &maxlen, Index nbElts, MemType memtype, Index &num_expansions)
Expand the existing storage.
RowVector3d w
Index memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t &glu)
Allocate various working space for the numerical factorization phase.
Point2 t(10, 10)
Index expand(VectorType &vec, Index &length, Index nbElts, Index keep_prev, Index &num_expansions)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:59