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>
|
---|
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>
|
---|
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>
|
---|
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>
|
---|
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+1) / n, m) * n; // estimated number of nonzeros in U
|
---|
157 | glu.nzlmax = (std::max)(Index(4), fillratio) * (annz+1) / 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>
|
---|
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
|
---|