Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
SparseLU_Memory.h
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
34namespace Eigen {
35namespace internal {
36
37enum { LUNoMarker = 3 };
38enum {emptyIdxLU = -1};
39inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b)
40{
41 return (std::max)(m, (t+b)*w);
42}
43
44template< typename Scalar>
45inline Index LUTempSpace(Index&m, Index& w)
46{
47 return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar);
48}
49
50
51
52
61template <typename Scalar, typename StorageIndex>
62template <typename VectorType>
63Index 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
150template <typename Scalar, typename StorageIndex>
151Index SparseLUImpl<Scalar,StorageIndex>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu)
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
207template <typename Scalar, typename StorageIndex>
208template <typename VectorType>
209Index 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
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74