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};
39 template<typename Index>
LUnumTempV(Index & m,Index & w,Index & t,Index & b)40 inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b)
41 {
42 return (std::max)(m, (t+b)*w);
43 }
44
45 template< typename Scalar, typename Index>
LUTempSpace(Index & m,Index & w)46 inline Index LUTempSpace(Index&m, Index& w)
47 {
48 return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar);
49 }
50
51
52
53
54 /**
55 * Expand the existing storage to accomodate more fill-ins
56 * \param vec Valid pointer to the vector to allocate or expand
57 * \param[in,out] length At input, contain the current length of the vector that is to be increased. At output, length of the newly allocated vector
58 * \param[in] nbElts Current number of elements in the factors
59 * \param keep_prev 1: use length and do not expand the vector; 0: compute new_len and expand
60 * \param[in,out] num_expansions Number of times the memory has been expanded
61 */
62 template <typename Scalar, typename Index>
63 template <typename VectorType>
expand(VectorType & vec,Index & length,Index nbElts,Index keep_prev,Index & num_expansions)64 Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions)
65 {
66
67 float alpha = 1.5; // Ratio of the memory increase
68 Index new_len; // New size of the allocated memory
69
70 if(num_expansions == 0 || keep_prev)
71 new_len = length ; // First time allocate requested
72 else
73 new_len = (std::max)(length+1,Index(alpha * length));
74
75 VectorType old_vec; // Temporary vector to hold the previous values
76 if (nbElts > 0 )
77 old_vec = vec.segment(0,nbElts);
78
79 //Allocate or expand the current vector
80 #ifdef EIGEN_EXCEPTIONS
81 try
82 #endif
83 {
84 vec.resize(new_len);
85 }
86 #ifdef EIGEN_EXCEPTIONS
87 catch(std::bad_alloc& )
88 #else
89 if(!vec.size())
90 #endif
91 {
92 if (!num_expansions)
93 {
94 // First time to allocate from LUMemInit()
95 // Let LUMemInit() deals with it.
96 return -1;
97 }
98 if (keep_prev)
99 {
100 // In this case, the memory length should not not be reduced
101 return new_len;
102 }
103 else
104 {
105 // Reduce the size and increase again
106 Index tries = 0; // Number of attempts
107 do
108 {
109 alpha = (alpha + 1)/2;
110 new_len = (std::max)(length+1,Index(alpha * length));
111 #ifdef EIGEN_EXCEPTIONS
112 try
113 #endif
114 {
115 vec.resize(new_len);
116 }
117 #ifdef EIGEN_EXCEPTIONS
118 catch(std::bad_alloc& )
119 #else
120 if (!vec.size())
121 #endif
122 {
123 tries += 1;
124 if ( tries > 10) return new_len;
125 }
126 } while (!vec.size());
127 }
128 }
129 //Copy the previous values to the newly allocated space
130 if (nbElts > 0)
131 vec.segment(0, nbElts) = old_vec;
132
133
134 length = new_len;
135 if(num_expansions) ++num_expansions;
136 return 0;
137 }
138
139 /**
140 * \brief Allocate various working space for the numerical factorization phase.
141 * \param m number of rows of the input matrix
142 * \param n number of columns
143 * \param annz number of initial nonzeros in the matrix
144 * \param lwork if lwork=-1, this routine returns an estimated size of the required memory
145 * \param glu persistent data to facilitate multiple factors : will be deleted later ??
146 * \param fillratio estimated ratio of fill in the factors
147 * \param panel_size Size of a panel
148 * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success
149 * \note Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation
150 */
151 template <typename Scalar, typename Index>
memInit(Index m,Index n,Index annz,Index lwork,Index fillratio,Index panel_size,GlobalLU_t & glu)152 Index SparseLUImpl<Scalar,Index>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu)
153 {
154 Index& num_expansions = glu.num_expansions; //No memory expansions so far
155 num_expansions = 0;
156 glu.nzumax = glu.nzlumax = (std::min)(fillratio * annz / n, m) * n; // estimated number of nonzeros in U
157 glu.nzlmax = (std::max)(Index(4), fillratio) * annz / 4; // estimated nnz in L factor
158 // Return the estimated size to the user if necessary
159 Index tempSpace;
160 tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar);
161 if (lwork == emptyIdxLU)
162 {
163 Index estimated_size;
164 estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace
165 + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) * sizeof(Scalar) + n;
166 return estimated_size;
167 }
168
169 // Setup the required space
170
171 // First allocate Integer pointers for L\U factors
172 glu.xsup.resize(n+1);
173 glu.supno.resize(n+1);
174 glu.xlsub.resize(n+1);
175 glu.xlusup.resize(n+1);
176 glu.xusub.resize(n+1);
177
178 // Reserve memory for L/U factors
179 do
180 {
181 if( (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0)
182 || (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0)
183 || (expand<IndexVector> (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0)
184 || (expand<IndexVector> (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) )
185 {
186 //Reduce the estimated size and retry
187 glu.nzlumax /= 2;
188 glu.nzumax /= 2;
189 glu.nzlmax /= 2;
190 if (glu.nzlumax < annz ) return glu.nzlumax;
191 }
192 } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
193
194 ++num_expansions;
195 return 0;
196
197 } // end LuMemInit
198
199 /**
200 * \brief Expand the existing storage
201 * \param vec vector to expand
202 * \param[in,out] maxlen On input, previous size of vec (Number of elements to copy ). on output, new size
203 * \param nbElts current number of elements in the vector.
204 * \param memtype Type of the element to expand
205 * \param num_expansions Number of expansions
206 * \return 0 on success, > 0 size of the memory allocated so far
207 */
208 template <typename Scalar, typename Index>
209 template <typename VectorType>
memXpand(VectorType & vec,Index & maxlen,Index nbElts,MemType memtype,Index & num_expansions)210 Index SparseLUImpl<Scalar,Index>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions)
211 {
212 Index failed_size;
213 if (memtype == USUB)
214 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions);
215 else
216 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions);
217
218 if (failed_size)
219 return failed_size;
220
221 return 0 ;
222 }
223
224 } // end namespace internal
225
226 } // end namespace Eigen
227 #endif // EIGEN_SPARSELU_MEMORY
228