|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === UMF_grow_front ======================================================= */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* -------------------------------------------------------------------------- */ |
|
|
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */ |
|
|
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */ |
|
|
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */ |
|
|
9 /* -------------------------------------------------------------------------- */ |
|
|
10 |
|
|
11 /* Current frontal matrix is too small. Make it bigger. */ |
|
|
12 |
|
|
13 #include "umf_internal.h" |
|
|
14 #include "umf_mem_free_tail_block.h" |
|
|
15 #include "umf_mem_alloc_tail_block.h" |
|
|
16 #include "umf_get_memory.h" |
|
|
17 |
|
|
18 GLOBAL Int UMF_grow_front |
|
|
19 ( |
|
|
20 NumericType *Numeric, |
|
|
21 Int fnr2, /* desired size is fnr2-by-fnc2 */ |
|
|
22 Int fnc2, |
|
|
23 WorkType *Work, |
|
|
24 Int do_what /* -1: UMF_start_front |
|
|
25 * 0: UMF_init_front, do not recompute Fcpos |
|
|
26 * 1: UMF_extend_front |
|
|
27 * 2: UMF_init_front, recompute Fcpos */ |
|
|
28 ) |
|
|
29 { |
|
|
30 /* ---------------------------------------------------------------------- */ |
|
|
31 /* local variables */ |
|
|
32 /* ---------------------------------------------------------------------- */ |
|
|
33 |
|
|
34 double s ; |
|
|
35 Entry *Fcold, *Fcnew ; |
|
|
36 Int j, i, col, *Fcpos, *Fcols, fnrows_max, fncols_max, fnr_curr, nb, |
|
|
37 fnrows_new, fncols_new, fnr_min, fnc_min, minsize, |
|
|
38 newsize, fnrows, fncols, *E, eloc ; |
|
|
39 |
|
|
40 /* ---------------------------------------------------------------------- */ |
|
|
41 /* get parameters */ |
|
|
42 /* ---------------------------------------------------------------------- */ |
|
|
43 |
|
|
44 #ifndef NDEBUG |
|
|
45 if (do_what != -1) UMF_debug++ ; |
|
|
46 DEBUG0 (("\n\n====================GROW FRONT: do_what: "ID"\n", do_what)) ; |
|
|
47 if (do_what != -1) UMF_debug-- ; |
|
|
48 ASSERT (Work->do_grow) ; |
|
|
49 ASSERT (Work->fnpiv == 0) ; |
|
|
50 #endif |
|
|
51 |
|
|
52 Fcols = Work->Fcols ; |
|
|
53 Fcpos = Work->Fcpos ; |
|
|
54 E = Work->E ; |
|
|
55 |
|
|
56 /* ---------------------------------------------------------------------- */ |
|
|
57 /* The current front is too small, find the new size */ |
|
|
58 /* ---------------------------------------------------------------------- */ |
|
|
59 |
|
|
60 /* maximum size of frontal matrix for this chain */ |
|
|
61 nb = Work->nb ; |
|
|
62 fnrows_max = Work->fnrows_max + nb ; |
|
|
63 fncols_max = Work->fncols_max + nb ; |
|
|
64 ASSERT (fnrows_max >= 0 && (fnrows_max % 2) == 1) ; |
|
|
65 DEBUG0 (("Max size: "ID"-by-"ID" (incl. "ID" pivot block\n", |
|
|
66 fnrows_max, fncols_max, nb)) ; |
|
|
67 |
|
|
68 /* current dimensions of frontal matrix: fnr-by-fnc */ |
|
|
69 DEBUG0 (("Current : "ID"-by-"ID" (excl "ID" pivot blocks)\n", |
|
|
70 Work->fnr_curr, Work->fnc_curr, nb)) ; |
|
|
71 ASSERT (Work->fnr_curr >= 0) ; |
|
|
72 ASSERT ((Work->fnr_curr % 2 == 1) || Work->fnr_curr == 0) ; |
|
|
73 |
|
|
74 /* required dimensions of frontal matrix: fnr_min-by-fnc_min */ |
|
|
75 fnrows_new = Work->fnrows_new + 1 ; |
|
|
76 fncols_new = Work->fncols_new + 1 ; |
|
|
77 ASSERT (fnrows_new >= 0) ; |
|
|
78 if (fnrows_new % 2 == 0) fnrows_new++ ; |
|
|
79 fnrows_new += nb ; |
|
|
80 fncols_new += nb ; |
|
|
81 fnr_min = MIN (fnrows_new, fnrows_max) ; |
|
|
82 fnc_min = MIN (fncols_new, fncols_max) ; |
|
|
83 minsize = fnr_min * fnc_min ; |
|
|
84 if (INT_OVERFLOW ((double) fnr_min * (double) fnc_min * sizeof (Entry))) |
|
|
85 { |
|
|
86 /* :: the minimum front size is bigger than the integer maximum :: */ |
|
|
87 return (FALSE) ; |
|
|
88 } |
|
|
89 ASSERT (fnr_min >= 0) ; |
|
|
90 ASSERT (fnr_min % 2 == 1) ; |
|
|
91 |
|
|
92 DEBUG0 (("Min : "ID"-by-"ID"\n", fnr_min, fnc_min)) ; |
|
|
93 |
|
|
94 /* grow the front to fnr2-by-fnc2, but no bigger than the maximum, |
|
|
95 * and no smaller than the minumum. */ |
|
|
96 DEBUG0 (("Desired : ("ID"+"ID")-by-("ID"+"ID")\n", fnr2, nb, fnc2, nb)) ; |
|
|
97 fnr2 += nb ; |
|
|
98 fnc2 += nb ; |
|
|
99 ASSERT (fnr2 >= 0) ; |
|
|
100 if (fnr2 % 2 == 0) fnr2++ ; |
|
|
101 fnr2 = MAX (fnr2, fnr_min) ; |
|
|
102 fnc2 = MAX (fnc2, fnc_min) ; |
|
|
103 fnr2 = MIN (fnr2, fnrows_max) ; |
|
|
104 fnc2 = MIN (fnc2, fncols_max) ; |
|
|
105 DEBUG0 (("Try : "ID"-by-"ID"\n", fnr2, fnc2)) ; |
|
|
106 ASSERT (fnr2 >= 0) ; |
|
|
107 ASSERT (fnr2 % 2 == 1) ; |
|
|
108 |
|
|
109 s = ((double) fnr2) * ((double) fnc2) ; |
|
|
110 if (INT_OVERFLOW (s * sizeof (Entry))) |
|
|
111 { |
|
|
112 /* :: frontal matrix size int overflow :: */ |
|
|
113 /* the desired front size is bigger than the integer maximum */ |
|
|
114 /* compute a such that a*a*s < Int_MAX / sizeof (Entry) */ |
|
|
115 double a = 0.9 * sqrt ((Int_MAX / sizeof (Entry)) / s) ; |
|
|
116 fnr2 = MAX (fnr_min, a * fnr2) ; |
|
|
117 fnc2 = MAX (fnc_min, a * fnc2) ; |
|
|
118 /* the new frontal size is a*r*a*c = a*a*s */ |
|
|
119 newsize = fnr2 * fnc2 ; |
|
|
120 ASSERT (fnr2 >= 0) ; |
|
|
121 if (fnr2 % 2 == 0) fnr2++ ; |
|
|
122 fnc2 = newsize / fnr2 ; |
|
|
123 } |
|
|
124 |
|
|
125 fnr2 = MAX (fnr2, fnr_min) ; |
|
|
126 fnc2 = MAX (fnc2, fnc_min) ; |
|
|
127 newsize = fnr2 * fnc2 ; |
|
|
128 |
|
|
129 ASSERT (fnr2 >= 0) ; |
|
|
130 ASSERT (fnr2 % 2 == 1) ; |
|
|
131 ASSERT (fnr2 >= fnr_min) ; |
|
|
132 ASSERT (fnc2 >= fnc_min) ; |
|
|
133 ASSERT (newsize >= minsize) ; |
|
|
134 |
|
|
135 /* ---------------------------------------------------------------------- */ |
|
|
136 /* free the current front if it is empty of any numerical values */ |
|
|
137 /* ---------------------------------------------------------------------- */ |
|
|
138 |
|
|
139 if (E [0] && do_what != 1) |
|
|
140 { |
|
|
141 /* free the current front, if it exists and has nothing in it */ |
|
|
142 DEBUG0 (("Freeing empty front\n")) ; |
|
|
143 UMF_mem_free_tail_block (Numeric, E [0]) ; |
|
|
144 E [0] = 0 ; |
|
|
145 Work->Flublock = (Entry *) NULL ; |
|
|
146 Work->Flblock = (Entry *) NULL ; |
|
|
147 Work->Fublock = (Entry *) NULL ; |
|
|
148 Work->Fcblock = (Entry *) NULL ; |
|
|
149 } |
|
|
150 |
|
|
151 /* ---------------------------------------------------------------------- */ |
|
|
152 /* allocate the new front, doing garbage collection if necessary */ |
|
|
153 /* ---------------------------------------------------------------------- */ |
|
|
154 |
|
|
155 #ifndef NDEBUG |
|
|
156 UMF_allocfail = FALSE ; |
|
|
157 if (UMF_gprob > 0) /* a double relop, but ignore NaN case */ |
|
|
158 { |
|
|
159 double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ; |
|
|
160 DEBUG1 (("Check random %e %e\n", rrr, UMF_gprob)) ; |
|
|
161 UMF_allocfail = rrr < UMF_gprob ; |
|
|
162 if (UMF_allocfail) DEBUGm2 (("Random garbage collection (grow)\n")) ; |
|
|
163 } |
|
|
164 #endif |
|
|
165 |
|
|
166 DEBUG0 (("Attempt size: "ID"-by-"ID"\n", fnr2, fnc2)) ; |
|
|
167 eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ; |
|
|
168 |
|
|
169 if (!eloc) |
|
|
170 { |
|
|
171 /* Do garbage collection, realloc, and try again. Compact the current |
|
|
172 * contribution block in the front to fnrows-by-fncols. Note that |
|
|
173 * there are no pivot rows/columns in current front. Do not recompute |
|
|
174 * Fcpos in UMF_garbage_collection. */ |
|
|
175 DEBUGm3 (("get_memory from umf_grow_front\n")) ; |
|
|
176 if (!UMF_get_memory (Numeric, Work, 1 + UNITS (Entry, newsize), |
|
|
177 Work->fnrows, Work->fncols, FALSE)) |
|
|
178 { |
|
|
179 /* :: out of memory in umf_grow_front :: */ |
|
|
180 return (FALSE) ; /* out of memory */ |
|
|
181 } |
|
|
182 DEBUG0 (("Attempt size: "ID"-by-"ID" again\n", fnr2, fnc2)) ; |
|
|
183 eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ; |
|
|
184 } |
|
|
185 |
|
|
186 /* try again with something smaller */ |
|
|
187 while ((fnr2 != fnr_min || fnc2 != fnc_min) && !eloc) |
|
|
188 { |
|
|
189 fnr2 = MIN (fnr2 - 2, fnr2 * UMF_REALLOC_REDUCTION) ; |
|
|
190 fnc2 = MIN (fnc2 - 2, fnc2 * UMF_REALLOC_REDUCTION) ; |
|
|
191 ASSERT (fnr_min >= 0) ; |
|
|
192 ASSERT (fnr_min % 2 == 1) ; |
|
|
193 fnr2 = MAX (fnr_min, fnr2) ; |
|
|
194 fnc2 = MAX (fnc_min, fnc2) ; |
|
|
195 ASSERT (fnr2 >= 0) ; |
|
|
196 if (fnr2 % 2 == 0) fnr2++ ; |
|
|
197 newsize = fnr2 * fnc2 ; |
|
|
198 DEBUGm3 (("Attempt smaller size: "ID"-by-"ID" minsize "ID"-by-"ID"\n", |
|
|
199 fnr2, fnc2, fnr_min, fnc_min)) ; |
|
|
200 eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ; |
|
|
201 } |
|
|
202 |
|
|
203 /* try again with the smallest possible size */ |
|
|
204 if (!eloc) |
|
|
205 { |
|
|
206 fnr2 = fnr_min ; |
|
|
207 fnc2 = fnc_min ; |
|
|
208 newsize = minsize ; |
|
|
209 DEBUG0 (("Attempt minsize: "ID"-by-"ID"\n", fnr2, fnc2)) ; |
|
|
210 eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ; |
|
|
211 } |
|
|
212 |
|
|
213 if (!eloc) |
|
|
214 { |
|
|
215 /* out of memory */ |
|
|
216 return (FALSE) ; |
|
|
217 } |
|
|
218 |
|
|
219 ASSERT (fnr2 >= 0) ; |
|
|
220 ASSERT (fnr2 % 2 == 1) ; |
|
|
221 ASSERT (fnr2 >= fnr_min && fnc2 >= fnc_min) ; |
|
|
222 |
|
|
223 /* ---------------------------------------------------------------------- */ |
|
|
224 /* copy the old frontal matrix into the new one */ |
|
|
225 /* ---------------------------------------------------------------------- */ |
|
|
226 |
|
|
227 /* old contribution block (if any) */ |
|
|
228 fnr_curr = Work->fnr_curr ; /* garbage collection can change fn*_curr */ |
|
|
229 ASSERT (fnr_curr >= 0) ; |
|
|
230 ASSERT ((fnr_curr % 2 == 1) || fnr_curr == 0) ; |
|
|
231 fnrows = Work->fnrows ; |
|
|
232 fncols = Work->fncols ; |
|
|
233 Fcold = Work->Fcblock ; |
|
|
234 |
|
|
235 /* remove nb from the sizes */ |
|
|
236 fnr2 -= nb ; |
|
|
237 fnc2 -= nb ; |
|
|
238 |
|
|
239 /* new frontal matrix */ |
|
|
240 Work->Flublock = (Entry *) (Numeric->Memory + eloc) ; |
|
|
241 Work->Flblock = Work->Flublock + nb * nb ; |
|
|
242 Work->Fublock = Work->Flblock + nb * fnr2 ; |
|
|
243 Work->Fcblock = Work->Fublock + nb * fnc2 ; |
|
|
244 Fcnew = Work->Fcblock ; |
|
|
245 |
|
|
246 if (E [0]) |
|
|
247 { |
|
|
248 /* copy the old contribution block into the new one */ |
|
|
249 for (j = 0 ; j < fncols ; j++) |
|
|
250 { |
|
|
251 col = Fcols [j] ; |
|
|
252 DEBUG1 (("copy col "ID" \n",col)) ; |
|
|
253 ASSERT (col >= 0 && col < Work->n_col) ; |
|
|
254 for (i = 0 ; i < fnrows ; i++) |
|
|
255 { |
|
|
256 Fcnew [i] = Fcold [i] ; |
|
|
257 } |
|
|
258 Fcnew += fnr2 ; |
|
|
259 Fcold += fnr_curr ; |
|
|
260 DEBUG1 (("new offset col "ID" "ID"\n",col, j * fnr2)) ; |
|
|
261 Fcpos [col] = j * fnr2 ; |
|
|
262 } |
|
|
263 } |
|
|
264 else if (do_what == 2) |
|
|
265 { |
|
|
266 /* just find the new column offsets */ |
|
|
267 for (j = 0 ; j < fncols ; j++) |
|
|
268 { |
|
|
269 col = Fcols [j] ; |
|
|
270 DEBUG1 (("new offset col "ID" "ID"\n",col, j * fnr2)) ; |
|
|
271 Fcpos [col] = j * fnr2 ; |
|
|
272 } |
|
|
273 } |
|
|
274 |
|
|
275 /* free the old frontal matrix */ |
|
|
276 UMF_mem_free_tail_block (Numeric, E [0]) ; |
|
|
277 |
|
|
278 /* ---------------------------------------------------------------------- */ |
|
|
279 /* new frontal matrix size */ |
|
|
280 /* ---------------------------------------------------------------------- */ |
|
|
281 |
|
|
282 E [0] = eloc ; |
|
|
283 Work->fnr_curr = fnr2 ; /* C block is fnr2-by-fnc2 */ |
|
|
284 Work->fnc_curr = fnc2 ; |
|
|
285 Work->fcurr_size = newsize ; /* including LU, L, U, and C blocks */ |
|
|
286 Work->do_grow = FALSE ; /* the front has just been grown */ |
|
|
287 |
|
|
288 ASSERT (Work->fnr_curr >= 0) ; |
|
|
289 ASSERT (Work->fnr_curr % 2 == 1) ; |
|
|
290 DEBUG0 (("Newly grown front: "ID"+"ID" by "ID"+"ID"\n", Work->fnr_curr, |
|
|
291 nb, Work->fnc_curr, nb)) ; |
|
|
292 return (TRUE) ; |
|
|
293 } |
