Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/Source/umf_grow_front.c @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
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 }