1 /* ========================================================================== */
2 /* === UMFPACK_numeric ====================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
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 /*
12 User-callable. Factorizes A into its LU factors, given a symbolic
13 pre-analysis computed by UMFPACK_symbolic. See umfpack_numeric.h for a
14 description.
15
16 Dynamic memory allocation: substantial. See comments (1) through (7),
17 below. If an error occurs, all allocated space is free'd by UMF_free.
18 If successful, the Numeric object contains 11 to 13 objects allocated by
19 UMF_malloc that hold the LU factors of the input matrix.
20 */
21
22 #include "umf_internal.h"
23 #include "umf_valid_symbolic.h"
24 #include "umf_set_stats.h"
25 #include "umf_kernel.h"
26 #include "umf_malloc.h"
27 #include "umf_free.h"
28 #include "umf_realloc.h"
29
30 #ifndef NDEBUG
31 PRIVATE Int init_count ;
32 #endif
33
34 PRIVATE Int work_alloc
35 (
36 WorkType *Work,
37 SymbolicType *Symbolic
38 ) ;
39
40 PRIVATE void free_work
41 (
42 WorkType *Work
43 ) ;
44
45 PRIVATE Int numeric_alloc
46 (
47 NumericType **NumericHandle,
48 SymbolicType *Symbolic,
49 double alloc_init,
50 Int scale
51 ) ;
52
53 PRIVATE void error
54 (
55 NumericType **Numeric,
56 WorkType *Work
57 ) ;
58
59
60 /* ========================================================================== */
61 /* === UMFPACK_numeric ====================================================== */
62 /* ========================================================================== */
63
UMFPACK_numeric(const Int Ap[],const Int Ai[],const double Ax[],const double Az[],void * SymbolicHandle,void ** NumericHandle,const double Control[UMFPACK_CONTROL],double User_Info[UMFPACK_INFO])64 GLOBAL Int UMFPACK_numeric
65 (
66 const Int Ap [ ],
67 const Int Ai [ ],
68 const double Ax [ ],
69 #ifdef COMPLEX
70 const double Az [ ],
71 #endif
72 void *SymbolicHandle,
73 void **NumericHandle,
74 const double Control [UMFPACK_CONTROL],
75 double User_Info [UMFPACK_INFO]
76 )
77 {
78
79 /* ---------------------------------------------------------------------- */
80 /* local variables */
81 /* ---------------------------------------------------------------------- */
82
83 double Info2 [UMFPACK_INFO], alloc_init, relpt, relpt2, droptol,
84 front_alloc_init, stats [2] ;
85 double *Info ;
86 WorkType WorkSpace, *Work ;
87 NumericType *Numeric ;
88 SymbolicType *Symbolic ;
89 Int n_row, n_col, n_inner, newsize, i, status, *inew, npiv, ulen, scale ;
90 Unit *mnew ;
91
92 /* ---------------------------------------------------------------------- */
93 /* get the amount of time used by the process so far */
94 /* ---------------------------------------------------------------------- */
95
96 umfpack_tic (stats) ;
97
98 /* ---------------------------------------------------------------------- */
99 /* initialize and check inputs */
100 /* ---------------------------------------------------------------------- */
101
102 #ifndef NDEBUG
103 UMF_dump_start ( ) ;
104 init_count = UMF_malloc_count ;
105 DEBUGm4 (("\nUMFPACK numeric: U transpose version\n")) ;
106 #endif
107
108 /* If front_alloc_init negative then allocate that size of front in
109 * UMF_start_front. If alloc_init negative, then allocate that initial
110 * size of Numeric->Memory. */
111
112 relpt = GET_CONTROL (UMFPACK_PIVOT_TOLERANCE,
113 UMFPACK_DEFAULT_PIVOT_TOLERANCE) ;
114 relpt2 = GET_CONTROL (UMFPACK_SYM_PIVOT_TOLERANCE,
115 UMFPACK_DEFAULT_SYM_PIVOT_TOLERANCE) ;
116 alloc_init = GET_CONTROL (UMFPACK_ALLOC_INIT, UMFPACK_DEFAULT_ALLOC_INIT) ;
117 front_alloc_init = GET_CONTROL (UMFPACK_FRONT_ALLOC_INIT,
118 UMFPACK_DEFAULT_FRONT_ALLOC_INIT) ;
119 scale = GET_CONTROL (UMFPACK_SCALE, UMFPACK_DEFAULT_SCALE) ;
120 droptol = GET_CONTROL (UMFPACK_DROPTOL, UMFPACK_DEFAULT_DROPTOL) ;
121
122 relpt = MAX (0.0, MIN (relpt, 1.0)) ;
123 relpt2 = MAX (0.0, MIN (relpt2, 1.0)) ;
124 droptol = MAX (0.0, droptol) ;
125 front_alloc_init = MIN (1.0, front_alloc_init) ;
126
127 if (scale != UMFPACK_SCALE_NONE && scale != UMFPACK_SCALE_MAX)
128 {
129 scale = UMFPACK_DEFAULT_SCALE ;
130 }
131
132 if (User_Info != (double *) NULL)
133 {
134 /* return Info in user's array */
135 Info = User_Info ;
136 /* clear the parts of Info that are set by UMFPACK_numeric */
137 for (i = UMFPACK_NUMERIC_SIZE ; i <= UMFPACK_MAX_FRONT_NCOLS ; i++)
138 {
139 Info [i] = EMPTY ;
140 }
141 for (i = UMFPACK_NUMERIC_DEFRAG ; i < UMFPACK_IR_TAKEN ; i++)
142 {
143 Info [i] = EMPTY ;
144 }
145 }
146 else
147 {
148 /* no Info array passed - use local one instead */
149 Info = Info2 ;
150 for (i = 0 ; i < UMFPACK_INFO ; i++)
151 {
152 Info [i] = EMPTY ;
153 }
154 }
155
156 Symbolic = (SymbolicType *) SymbolicHandle ;
157 Numeric = (NumericType *) NULL ;
158 if (!UMF_valid_symbolic (Symbolic))
159 {
160 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_Symbolic_object ;
161 return (UMFPACK_ERROR_invalid_Symbolic_object) ;
162 }
163
164 /* compute alloc_init automatically for AMD ordering */
165 if (Symbolic->ordering == UMFPACK_ORDERING_AMD && alloc_init >= 0)
166 {
167 alloc_init = (Symbolic->nz + Symbolic->amd_lunz) / Symbolic->lunz_bound;
168 alloc_init = MIN (1.0, alloc_init) ;
169 alloc_init *= UMF_REALLOC_INCREASE ;
170 }
171
172 n_row = Symbolic->n_row ;
173 n_col = Symbolic->n_col ;
174 n_inner = MIN (n_row, n_col) ;
175
176 /* check for integer overflow in Numeric->Memory minimum size */
177 if (INT_OVERFLOW (Symbolic->dnum_mem_init_usage * sizeof (Unit)))
178 {
179 /* :: int overflow, initial Numeric->Memory size :: */
180 /* There's no hope to allocate a Numeric object big enough simply to
181 * hold the initial matrix, so return an out-of-memory condition */
182 DEBUGm4 (("out of memory: numeric int overflow\n")) ;
183 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
184 return (UMFPACK_ERROR_out_of_memory) ;
185 }
186
187 Info [UMFPACK_STATUS] = UMFPACK_OK ;
188 Info [UMFPACK_NROW] = n_row ;
189 Info [UMFPACK_NCOL] = n_col ;
190 Info [UMFPACK_SIZE_OF_UNIT] = (double) (sizeof (Unit)) ;
191
192 if (!Ap || !Ai || !Ax || !NumericHandle)
193 {
194 Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
195 return (UMFPACK_ERROR_argument_missing) ;
196 }
197
198 Info [UMFPACK_NZ] = Ap [n_col] ;
199 *NumericHandle = (void *) NULL ;
200
201 /* ---------------------------------------------------------------------- */
202 /* allocate the Work object */
203 /* ---------------------------------------------------------------------- */
204
205 /* (1) calls UMF_malloc 15 or 17 times, to obtain temporary workspace of
206 * size c+1 Entry's and 2*(n_row+1) + 3*(n_col+1) + (n_col+n_inner+1) +
207 * (nn+1) + * 3*(c+1) + 2*(r+1) + max(r,c) + (nfr+1) integers plus 2*nn
208 * more integers if diagonal pivoting is to be done. r is the maximum
209 * number of rows in any frontal matrix, c is the maximum number of columns
210 * in any frontal matrix, n_inner is min (n_row,n_col), nn is
211 * max (n_row,n_col), and nfr is the number of frontal matrices. For a
212 * square matrix, this is c+1 Entry's and about 8n + 3c + 2r + max(r,c) +
213 * nfr integers, plus 2n more for diagonal pivoting.
214 */
215
216 Work = &WorkSpace ;
217 Work->n_row = n_row ;
218 Work->n_col = n_col ;
219 Work->nfr = Symbolic->nfr ;
220 Work->nb = Symbolic->nb ;
221 Work->n1 = Symbolic->n1 ;
222
223 if (!work_alloc (Work, Symbolic))
224 {
225 DEBUGm4 (("out of memory: numeric work\n")) ;
226 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
227 error (&Numeric, Work) ;
228 return (UMFPACK_ERROR_out_of_memory) ;
229 }
230 ASSERT (UMF_malloc_count == init_count + 16 + 2*Symbolic->prefer_diagonal) ;
231
232 /* ---------------------------------------------------------------------- */
233 /* allocate Numeric object */
234 /* ---------------------------------------------------------------------- */
235
236 /* (2) calls UMF_malloc 10 or 11 times, for a total space of
237 * sizeof (NumericType) bytes, 4*(n_row+1) + 4*(n_row+1) integers, and
238 * (n_inner+1) Entry's, plus n_row Entry's if row scaling is to be done.
239 * sizeof (NumericType) is a small constant. Next, it calls UMF_malloc
240 * once, for the variable-sized part of the Numeric object
241 * (Numeric->Memory). The size of this object is the larger of
242 * (Control [UMFPACK_ALLOC_INIT]) * (the approximate upper bound computed
243 * by UMFPACK_symbolic), and the minimum required to start the numerical
244 * factorization. * This request is reduced if it fails.
245 */
246
247 if (!numeric_alloc (&Numeric, Symbolic, alloc_init, scale))
248 {
249 DEBUGm4 (("out of memory: initial numeric\n")) ;
250 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
251 error (&Numeric, Work) ;
252 return (UMFPACK_ERROR_out_of_memory) ;
253 }
254 DEBUG0 (("malloc: init_count " ID " UMF_malloc_count " ID "\n",
255 init_count, UMF_malloc_count)) ;
256 ASSERT (UMF_malloc_count == init_count
257 + (16 + 2*Symbolic->prefer_diagonal)
258 + (11 + (scale != UMFPACK_SCALE_NONE))) ;
259
260 /* set control parameters */
261 Numeric->relpt = relpt ;
262 Numeric->relpt2 = relpt2 ;
263 Numeric->droptol = droptol ;
264 Numeric->alloc_init = alloc_init ;
265 Numeric->front_alloc_init = front_alloc_init ;
266 Numeric->scale = scale ;
267
268 DEBUG0 (("umf relpt %g %g init %g %g inc %g red %g\n",
269 relpt, relpt2, alloc_init, front_alloc_init,
270 UMF_REALLOC_INCREASE, UMF_REALLOC_REDUCTION)) ;
271
272 /* ---------------------------------------------------------------------- */
273 /* scale and factorize */
274 /* ---------------------------------------------------------------------- */
275
276 /* (3) During numerical factorization (inside UMF_kernel), the variable-size
277 * block of memory is increased in size via a call to UMF_realloc if it is
278 * found to be too small. During factorization, this block holds the
279 * pattern and values of L and U at the top end, and the elements
280 * (contibution blocks) and the current frontal matrix (Work->F*) at the
281 * bottom end. The peak size of the variable-sized object is estimated in
282 * UMFPACK_*symbolic (Info [UMFPACK_VARIABLE_PEAK_ESTIMATE]), although this
283 * upper bound can be very loose. The size of the Symbolic object
284 * (which is currently allocated) is in Info [UMFPACK_SYMBOLIC_SIZE], and
285 * is between 2*n and 13*n integers.
286 */
287
288 DEBUG0 (("Calling umf_kernel\n")) ;
289 status = UMF_kernel (Ap, Ai, Ax,
290 #ifdef COMPLEX
291 Az,
292 #endif
293 Numeric, Work, Symbolic) ;
294
295 Info [UMFPACK_STATUS] = status ;
296 if (status < UMFPACK_OK)
297 {
298 /* out of memory, or pattern has changed */
299 error (&Numeric, Work) ;
300 return (status) ;
301 }
302
303 Info [UMFPACK_FORCED_UPDATES] = Work->nforced ;
304 Info [UMFPACK_VARIABLE_INIT] = Numeric->init_usage ;
305 if (Symbolic->prefer_diagonal)
306 {
307 Info [UMFPACK_NOFF_DIAG] = Work->noff_diagonal ;
308 }
309
310 DEBUG0 (("malloc: init_count " ID " UMF_malloc_count " ID "\n",
311 init_count, UMF_malloc_count)) ;
312
313 npiv = Numeric->npiv ; /* = n_inner for nonsingular matrices */
314 ulen = Numeric->ulen ; /* = 0 for square nonsingular matrices */
315
316 /* ---------------------------------------------------------------------- */
317 /* free Work object */
318 /* ---------------------------------------------------------------------- */
319
320 /* (4) After numerical factorization all of the objects allocated in step
321 * (1) are freed via UMF_free, except that one object of size n_col+1 is
322 * kept if there are off-diagonal nonzeros in the last pivot row (can only
323 * occur for singular or rectangular matrices). This is Work->Upattern,
324 * which is transfered to Numeric->Upattern if ulen > 0.
325 */
326
327 DEBUG0 (("malloc: init_count " ID " UMF_malloc_count " ID "\n",
328 init_count, UMF_malloc_count)) ;
329
330 free_work (Work) ;
331
332 DEBUG0 (("malloc: init_count " ID " UMF_malloc_count " ID "\n",
333 init_count, UMF_malloc_count)) ;
334 DEBUG0 (("Numeric->ulen: " ID " scale: " ID "\n", ulen, scale)) ;
335 ASSERT (UMF_malloc_count == init_count + (ulen > 0) +
336 (11 + (scale != UMFPACK_SCALE_NONE))) ;
337
338 /* ---------------------------------------------------------------------- */
339 /* reduce Lpos, Lilen, Lip, Upos, Uilen and Uip to size npiv+1 */
340 /* ---------------------------------------------------------------------- */
341
342 /* (5) Six components of the Numeric object are reduced in size if the
343 * matrix is singular or rectangular. The original size is 3*(n_row+1) +
344 * 3*(n_col+1) integers. The new size is 6*(npiv+1) integers. For
345 * square non-singular matrices, these two sizes are the same.
346 */
347
348 if (npiv < n_row)
349 {
350 /* reduce Lpos, Uilen, and Uip from size n_row+1 to size npiv */
351 inew = (Int *) UMF_realloc (Numeric->Lpos, npiv+1, sizeof (Int)) ;
352 if (inew)
353 {
354 Numeric->Lpos = inew ;
355 }
356 inew = (Int *) UMF_realloc (Numeric->Uilen, npiv+1, sizeof (Int)) ;
357 if (inew)
358 {
359 Numeric->Uilen = inew ;
360 }
361 inew = (Int *) UMF_realloc (Numeric->Uip, npiv+1, sizeof (Int)) ;
362 if (inew)
363 {
364 Numeric->Uip = inew ;
365 }
366 }
367
368 if (npiv < n_col)
369 {
370 /* reduce Upos, Lilen, and Lip from size n_col+1 to size npiv */
371 inew = (Int *) UMF_realloc (Numeric->Upos, npiv+1, sizeof (Int)) ;
372 if (inew)
373 {
374 Numeric->Upos = inew ;
375 }
376 inew = (Int *) UMF_realloc (Numeric->Lilen, npiv+1, sizeof (Int)) ;
377 if (inew)
378 {
379 Numeric->Lilen = inew ;
380 }
381 inew = (Int *) UMF_realloc (Numeric->Lip, npiv+1, sizeof (Int)) ;
382 if (inew)
383 {
384 Numeric->Lip = inew ;
385 }
386 }
387
388 /* ---------------------------------------------------------------------- */
389 /* reduce Numeric->Upattern from size n_col+1 to size ulen+1 */
390 /* ---------------------------------------------------------------------- */
391
392 /* (6) The size of Numeric->Upattern (formerly Work->Upattern) is reduced
393 * from size n_col+1 to size ulen + 1. If ulen is zero, the object does
394 * not exist. */
395
396 DEBUG4 (("ulen: " ID " Upattern " ID "\n", ulen, (Int) Numeric->Upattern)) ;
397 ASSERT (IMPLIES (ulen == 0, Numeric->Upattern == (Int *) NULL)) ;
398 if (ulen > 0 && ulen < n_col)
399 {
400 inew = (Int *) UMF_realloc (Numeric->Upattern, ulen+1, sizeof (Int)) ;
401 if (inew)
402 {
403 Numeric->Upattern = inew ;
404 }
405 }
406
407 /* ---------------------------------------------------------------------- */
408 /* reduce Numeric->Memory to hold just the LU factors at the head */
409 /* ---------------------------------------------------------------------- */
410
411 /* (7) The variable-sized block (Numeric->Memory) is reduced to hold just L
412 * and U, via a call to UMF_realloc, since the frontal matrices are no
413 * longer needed.
414 */
415
416 newsize = Numeric->ihead ;
417 if (newsize < Numeric->size)
418 {
419 mnew = (Unit *) UMF_realloc (Numeric->Memory, newsize, sizeof (Unit)) ;
420 if (mnew)
421 {
422 /* realloc succeeded (how can it fail since the size is reduced?) */
423 Numeric->Memory = mnew ;
424 Numeric->size = newsize ;
425 }
426 }
427 Numeric->ihead = Numeric->size ;
428 Numeric->itail = Numeric->ihead ;
429 Numeric->tail_usage = 0 ;
430 Numeric->ibig = EMPTY ;
431 /* UMF_mem_alloc_tail_block can no longer be called (no tail marker) */
432
433 /* ---------------------------------------------------------------------- */
434 /* report the results and return the Numeric object */
435 /* ---------------------------------------------------------------------- */
436
437 UMF_set_stats (
438 Info,
439 Symbolic,
440 (double) Numeric->max_usage, /* actual peak Numeric->Memory */
441 (double) Numeric->size, /* actual final Numeric->Memory */
442 Numeric->flops, /* actual "true flops" */
443 (double) Numeric->lnz + n_inner, /* actual nz in L */
444 (double) Numeric->unz + Numeric->nnzpiv, /* actual nz in U */
445 (double) Numeric->maxfrsize, /* actual largest front size */
446 (double) ulen, /* actual Numeric->Upattern size */
447 (double) npiv, /* actual # pivots found */
448 (double) Numeric->maxnrows, /* actual largest #rows in front */
449 (double) Numeric->maxncols, /* actual largest #cols in front */
450 scale != UMFPACK_SCALE_NONE,
451 Symbolic->prefer_diagonal,
452 ACTUAL) ;
453
454 Info [UMFPACK_ALLOC_INIT_USED] = Numeric->alloc_init ;
455 Info [UMFPACK_NUMERIC_DEFRAG] = Numeric->ngarbage ;
456 Info [UMFPACK_NUMERIC_REALLOC] = Numeric->nrealloc ;
457 Info [UMFPACK_NUMERIC_COSTLY_REALLOC] = Numeric->ncostly ;
458 Info [UMFPACK_COMPRESSED_PATTERN] = Numeric->isize ;
459 Info [UMFPACK_LU_ENTRIES] = Numeric->nLentries + Numeric->nUentries +
460 Numeric->npiv ;
461 Info [UMFPACK_UDIAG_NZ] = Numeric->nnzpiv ;
462 Info [UMFPACK_RSMIN] = Numeric->rsmin ;
463 Info [UMFPACK_RSMAX] = Numeric->rsmax ;
464 Info [UMFPACK_WAS_SCALED] = Numeric->scale ;
465
466 /* nz in L and U with no dropping of small entries */
467 Info [UMFPACK_ALL_LNZ] = Numeric->all_lnz + n_inner ;
468 Info [UMFPACK_ALL_UNZ] = Numeric->all_unz + Numeric->nnzpiv ;
469 Info [UMFPACK_NZDROPPED] =
470 (Numeric->all_lnz - Numeric->lnz)
471 + (Numeric->all_unz - Numeric->unz) ;
472
473 /* estimate of the reciprocal of the condition number. */
474 if (SCALAR_IS_ZERO (Numeric->min_udiag)
475 || SCALAR_IS_ZERO (Numeric->max_udiag)
476 || SCALAR_IS_NAN (Numeric->min_udiag)
477 || SCALAR_IS_NAN (Numeric->max_udiag))
478 {
479 /* rcond is zero if there is any zero or NaN on the diagonal */
480 Numeric->rcond = 0.0 ;
481 }
482 else
483 {
484 /* estimate of the recipricol of the condition number. */
485 /* This is NaN if diagonal is zero-free, but has one or more NaN's. */
486 Numeric->rcond = Numeric->min_udiag / Numeric->max_udiag ;
487 }
488 Info [UMFPACK_UMIN] = Numeric->min_udiag ;
489 Info [UMFPACK_UMAX] = Numeric->max_udiag ;
490 Info [UMFPACK_RCOND] = Numeric->rcond ;
491
492 if (Numeric->nnzpiv < n_inner
493 || SCALAR_IS_ZERO (Numeric->rcond) || SCALAR_IS_NAN (Numeric->rcond))
494 {
495 /* there are zeros and/or NaN's on the diagonal of U */
496 DEBUG0 (("Warning, matrix is singular in umfpack_numeric\n")) ;
497 DEBUG0 (("nnzpiv " ID " n_inner " ID " rcond %g\n", Numeric->nnzpiv,
498 n_inner, Numeric->rcond)) ;
499 status = UMFPACK_WARNING_singular_matrix ;
500 Info [UMFPACK_STATUS] = status ;
501 }
502
503 Numeric->valid = NUMERIC_VALID ;
504 *NumericHandle = (void *) Numeric ;
505
506 /* Numeric has 11 to 13 objects */
507 ASSERT (UMF_malloc_count == init_count + 11 +
508 + (ulen > 0) /* Numeric->Upattern */
509 + (scale != UMFPACK_SCALE_NONE)) ; /* Numeric->Rs */
510
511 /* ---------------------------------------------------------------------- */
512 /* get the time used by UMFPACK_numeric */
513 /* ---------------------------------------------------------------------- */
514
515 umfpack_toc (stats) ;
516 Info [UMFPACK_NUMERIC_WALLTIME] = stats [0] ;
517 Info [UMFPACK_NUMERIC_TIME] = stats [1] ;
518
519 /* return UMFPACK_OK or UMFPACK_WARNING_singular_matrix */
520 return (status) ;
521
522 }
523
524
525 /* ========================================================================== */
526 /* === numeric_alloc ======================================================== */
527 /* ========================================================================== */
528
529 /* Allocate the Numeric object */
530
numeric_alloc(NumericType ** NumericHandle,SymbolicType * Symbolic,double alloc_init,Int scale)531 PRIVATE Int numeric_alloc
532 (
533 NumericType **NumericHandle,
534 SymbolicType *Symbolic,
535 double alloc_init,
536 Int scale
537 )
538 {
539 double nsize, bsize ;
540 Int n_row, n_col, n_inner, min_usage, trying ;
541 NumericType *Numeric ;
542
543 DEBUG0 (("numeric alloc:\n")) ;
544
545 n_row = Symbolic->n_row ;
546 n_col = Symbolic->n_col ;
547 n_inner = MIN (n_row, n_col) ;
548 *NumericHandle = (NumericType *) NULL ;
549
550 /* 1 allocation: accounted for in UMF_set_stats (num_On_size1),
551 * free'd in umfpack_free_numeric */
552 Numeric = (NumericType *) UMF_malloc (1, sizeof (NumericType)) ;
553
554 if (!Numeric)
555 {
556 return (FALSE) ; /* out of memory */
557 }
558 Numeric->valid = 0 ;
559 *NumericHandle = Numeric ;
560
561 /* 9 allocations: accounted for in UMF_set_stats (num_On_size1),
562 * free'd in umfpack_free_numeric */
563 Numeric->D = (Entry *) UMF_malloc (n_inner+1, sizeof (Entry)) ;
564 Numeric->Rperm = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
565 Numeric->Cperm = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
566 Numeric->Lpos = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
567 Numeric->Lilen = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
568 Numeric->Lip = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
569 Numeric->Upos = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
570 Numeric->Uilen = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
571 Numeric->Uip = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
572
573 /* 1 allocation if scaling: in UMF_set_stats (num_On_size1),
574 * free'd in umfpack_free_numeric */
575 if (scale != UMFPACK_SCALE_NONE)
576 {
577 DEBUG0 (("Allocating scale factors\n")) ;
578 Numeric->Rs = (double *) UMF_malloc (n_row, sizeof (double)) ;
579 }
580 else
581 {
582 DEBUG0 (("No scale factors allocated (R = I)\n")) ;
583 Numeric->Rs = (double *) NULL ;
584 }
585
586 Numeric->Memory = (Unit *) NULL ;
587
588 /* Upattern has already been allocated as part of the Work object. If
589 * the matrix is singular or rectangular, and there are off-diagonal
590 * nonzeros in the last pivot row, then Work->Upattern is not free'd.
591 * Instead it is transfered to Numeric->Upattern. If it exists,
592 * Numeric->Upattern is free'd in umfpack_free_numeric. */
593 Numeric->Upattern = (Int *) NULL ; /* used for singular matrices only */
594
595 if (!Numeric->D || !Numeric->Rperm || !Numeric->Cperm || !Numeric->Upos ||
596 !Numeric->Lpos || !Numeric->Lilen || !Numeric->Uilen || !Numeric->Lip ||
597 !Numeric->Uip || (scale != UMFPACK_SCALE_NONE && !Numeric->Rs))
598 {
599 return (FALSE) ; /* out of memory */
600 }
601
602 /* ---------------------------------------------------------------------- */
603 /* allocate initial Numeric->Memory for LU factors and elements */
604 /* ---------------------------------------------------------------------- */
605
606 if (alloc_init < 0)
607 {
608 /* -alloc_init is the exact size to initially allocate */
609 nsize = -alloc_init ;
610 }
611 else
612 {
613 /* alloc_init is a ratio of the upper bound memory usage */
614 nsize = (alloc_init * Symbolic->num_mem_usage_est) + 1 ;
615 }
616 min_usage = Symbolic->num_mem_init_usage ;
617
618 /* Numeric->Memory must be large enough for UMF_kernel_init */
619 nsize = MAX (min_usage, nsize) ;
620
621 /* Numeric->Memory cannot be larger in size than Int_MAX / sizeof(Unit) */
622 /* For ILP32 mode: 2GB (nsize cannot be bigger than 256 Mwords) */
623 bsize = ((double) Int_MAX) / sizeof (Unit) - 1 ;
624 DEBUG0 (("bsize %g\n", bsize)) ;
625 nsize = MIN (nsize, bsize) ;
626
627 Numeric->size = (Int) nsize ;
628
629 DEBUG0 (("Num init %g usage_est %g numsize " ID " minusage " ID "\n",
630 alloc_init, Symbolic->num_mem_usage_est, Numeric->size, min_usage)) ;
631
632 /* allocates 1 object: */
633 /* keep trying until successful, or memory request is too small */
634 trying = TRUE ;
635 while (trying)
636 {
637 Numeric->Memory = (Unit *) UMF_malloc (Numeric->size, sizeof (Unit)) ;
638 if (Numeric->Memory)
639 {
640 DEBUG0 (("Successful Numeric->size: " ID "\n", Numeric->size)) ;
641 return (TRUE) ;
642 }
643 /* too much, reduce the request (but not below the minimum) */
644 /* and try again */
645 trying = Numeric->size > min_usage ;
646 Numeric->size = (Int)
647 (UMF_REALLOC_REDUCTION * ((double) Numeric->size)) ;
648 Numeric->size = MAX (min_usage, Numeric->size) ;
649 }
650
651 return (FALSE) ; /* we failed to allocate Numeric->Memory */
652 }
653
654
655 /* ========================================================================== */
656 /* === work_alloc =========================================================== */
657 /* ========================================================================== */
658
659 /* Allocate the Work object. Return TRUE if successful. */
660
work_alloc(WorkType * Work,SymbolicType * Symbolic)661 PRIVATE Int work_alloc
662 (
663 WorkType *Work,
664 SymbolicType *Symbolic
665 )
666 {
667 Int n_row, n_col, nn, maxnrows, maxncols, nfr, ok, maxnrc, n1 ;
668
669 n_row = Work->n_row ;
670 n_col = Work->n_col ;
671 nn = MAX (n_row, n_col) ;
672 nfr = Work->nfr ;
673 n1 = Symbolic->n1 ;
674 ASSERT (n1 <= n_row && n1 <= n_col) ;
675
676 maxnrows = Symbolic->maxnrows + Symbolic->nb ;
677 maxnrows = MIN (n_row, maxnrows) ;
678 maxncols = Symbolic->maxncols + Symbolic->nb ;
679 maxncols = MIN (n_col, maxncols) ;
680 maxnrc = MAX (maxnrows, maxncols) ;
681
682 DEBUG0 (("work alloc: maxnrows+nb " ID " maxncols+nb " ID "\n",
683 maxnrows, maxncols)) ;
684
685 /* 15 allocations, freed in free_work: */
686 /* accounted for in UMF_set_stats (work_usage) */
687 Work->Wx = (Entry *) UMF_malloc (maxnrows + 1, sizeof (Entry)) ;
688 Work->Wy = (Entry *) UMF_malloc (maxnrows + 1, sizeof (Entry)) ;
689 Work->Frpos = (Int *) UMF_malloc (n_row + 1, sizeof (Int)) ;
690 Work->Lpattern = (Int *) UMF_malloc (n_row + 1, sizeof (Int)) ;
691 Work->Fcpos = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
692 Work->Wp = (Int *) UMF_malloc (nn + 1, sizeof (Int)) ;
693 Work->Wrp = (Int *) UMF_malloc (MAX (n_col,maxnrows) + 1, sizeof (Int)) ;
694 Work->Frows = (Int *) UMF_malloc (maxnrows + 1, sizeof (Int)) ;
695 Work->Wm = (Int *) UMF_malloc (maxnrows + 1, sizeof (Int)) ;
696 Work->Fcols = (Int *) UMF_malloc (maxncols + 1, sizeof (Int)) ;
697 Work->Wio = (Int *) UMF_malloc (maxncols + 1, sizeof (Int)) ;
698 Work->Woi = (Int *) UMF_malloc (maxncols + 1, sizeof (Int)) ;
699 Work->Woo = (Int *) UMF_malloc (maxnrc + 1, sizeof (Int));
700 Work->elen = (n_col - n1) + (n_row - n1) + MIN (n_col-n1, n_row-n1) + 1 ;
701 Work->E = (Int *) UMF_malloc (Work->elen, sizeof (Int)) ;
702 Work->Front_new1strow = (Int *) UMF_malloc (nfr + 1, sizeof (Int)) ;
703
704 ok = (Work->Frpos && Work->Fcpos && Work->Lpattern
705 && Work->Wp && Work->Wrp && Work->Frows && Work->Fcols
706 && Work->Wio && Work->Woi && Work->Woo && Work->Wm
707 && Work->E && Work->Front_new1strow && Work->Wx && Work->Wy) ;
708
709 /* 2 allocations: accounted for in UMF_set_stats (work_usage) */
710 if (Symbolic->prefer_diagonal)
711 {
712 Work->Diagonal_map = (Int *) UMF_malloc (nn, sizeof (Int)) ;
713 Work->Diagonal_imap = (Int *) UMF_malloc (nn, sizeof (Int)) ;
714 ok = ok && Work->Diagonal_map && Work->Diagonal_imap ;
715 }
716 else
717 {
718 /* no diagonal map needed for rectangular matrices */
719 Work->Diagonal_map = (Int *) NULL ;
720 Work->Diagonal_imap = (Int *) NULL ;
721 }
722
723 /* 1 allocation, may become part of Numeric (if singular or rectangular): */
724 Work->Upattern = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
725 ok = ok && Work->Upattern ;
726
727 /* current frontal matrix does not yet exist */
728 Work->Flublock = (Entry *) NULL ;
729 Work->Flblock = (Entry *) NULL ;
730 Work->Fublock = (Entry *) NULL ;
731 Work->Fcblock = (Entry *) NULL ;
732
733 DEBUG0 (("work alloc done.\n")) ;
734 return (ok) ;
735 }
736
737
738 /* ========================================================================== */
739 /* === free_work ============================================================ */
740 /* ========================================================================== */
741
free_work(WorkType * Work)742 PRIVATE void free_work
743 (
744 WorkType *Work
745 )
746 {
747 DEBUG0 (("work free:\n")) ;
748 if (Work)
749 {
750 /* these 16 objects do exist */
751 Work->Wx = (Entry *) UMF_free ((void *) Work->Wx) ;
752 Work->Wy = (Entry *) UMF_free ((void *) Work->Wy) ;
753 Work->Frpos = (Int *) UMF_free ((void *) Work->Frpos) ;
754 Work->Fcpos = (Int *) UMF_free ((void *) Work->Fcpos) ;
755 Work->Lpattern = (Int *) UMF_free ((void *) Work->Lpattern) ;
756 Work->Upattern = (Int *) UMF_free ((void *) Work->Upattern) ;
757 Work->Wp = (Int *) UMF_free ((void *) Work->Wp) ;
758 Work->Wrp = (Int *) UMF_free ((void *) Work->Wrp) ;
759 Work->Frows = (Int *) UMF_free ((void *) Work->Frows) ;
760 Work->Fcols = (Int *) UMF_free ((void *) Work->Fcols) ;
761 Work->Wio = (Int *) UMF_free ((void *) Work->Wio) ;
762 Work->Woi = (Int *) UMF_free ((void *) Work->Woi) ;
763 Work->Woo = (Int *) UMF_free ((void *) Work->Woo) ;
764 Work->Wm = (Int *) UMF_free ((void *) Work->Wm) ;
765 Work->E = (Int *) UMF_free ((void *) Work->E) ;
766 Work->Front_new1strow =
767 (Int *) UMF_free ((void *) Work->Front_new1strow) ;
768
769 /* these objects might not exist */
770 Work->Diagonal_map = (Int *) UMF_free ((void *) Work->Diagonal_map) ;
771 Work->Diagonal_imap = (Int *) UMF_free ((void *) Work->Diagonal_imap) ;
772 }
773 DEBUG0 (("work free done.\n")) ;
774 }
775
776
777 /* ========================================================================== */
778 /* === error ================================================================ */
779 /* ========================================================================== */
780
781 /* Error return from UMFPACK_numeric. Free all allocated memory. */
782
error(NumericType ** Numeric,WorkType * Work)783 PRIVATE void error
784 (
785 NumericType **Numeric,
786 WorkType *Work
787 )
788 {
789 free_work (Work) ;
790 UMFPACK_free_numeric ((void **) Numeric) ;
791 ASSERT (UMF_malloc_count == init_count) ;
792 }
793