1 /* ========================================================================== */
2 /* === UMFPACK_numeric ====================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /*
11 User-callable. Factorizes A into its LU factors, given a symbolic
12 pre-analysis computed by UMFPACK_symbolic. See umfpack_numeric.h for a
13 description.
14
15 Dynamic memory allocation: substantial. See comments (1) through (7),
16 below. If an error occurs, all allocated space is free'd by UMF_free.
17 If successful, the Numeric object contains 11 to 13 objects allocated by
18 UMF_malloc that hold the LU factors of the input matrix.
19 */
20
21 #include "umf_internal.h"
22 #include "umf_valid_symbolic.h"
23 #include "umf_set_stats.h"
24 #include "umf_kernel.h"
25 #include "umf_malloc.h"
26 #include "umf_free.h"
27 #include "umf_realloc.h"
28
29 #ifndef NDEBUG
30 PRIVATE Int init_count ;
31 #endif
32
33 PRIVATE Int work_alloc
34 (
35 WorkType *Work,
36 SymbolicType *Symbolic
37 ) ;
38
39 PRIVATE void free_work
40 (
41 WorkType *Work
42 ) ;
43
44 PRIVATE Int numeric_alloc
45 (
46 NumericType **NumericHandle,
47 SymbolicType *Symbolic,
48 double alloc_init,
49 Int scale
50 ) ;
51
52 PRIVATE void error
53 (
54 NumericType **Numeric,
55 WorkType *Work
56 ) ;
57
58
59 /* ========================================================================== */
60 /* === UMFPACK_numeric ====================================================== */
61 /* ========================================================================== */
62
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])63 GLOBAL Int UMFPACK_numeric
64 (
65 const Int Ap [ ],
66 const Int Ai [ ],
67 const double Ax [ ],
68 #ifdef COMPLEX
69 const double Az [ ],
70 #endif
71 void *SymbolicHandle,
72 void **NumericHandle,
73 const double Control [UMFPACK_CONTROL],
74 double User_Info [UMFPACK_INFO]
75 )
76 {
77
78 /* ---------------------------------------------------------------------- */
79 /* local variables */
80 /* ---------------------------------------------------------------------- */
81
82 double Info2 [UMFPACK_INFO], alloc_init, relpt, relpt2, droptol,
83 front_alloc_init, stats [2] ;
84 double *Info ;
85 WorkType WorkSpace, *Work ;
86 NumericType *Numeric ;
87 SymbolicType *Symbolic ;
88 Int n_row, n_col, n_inner, newsize, i, status, *inew, npiv, ulen, scale ;
89 Unit *mnew ;
90
91 /* ---------------------------------------------------------------------- */
92 /* get the amount of time used by the process so far */
93 /* ---------------------------------------------------------------------- */
94
95 umfpack_tic (stats) ;
96
97 /* ---------------------------------------------------------------------- */
98 /* initialize and check inputs */
99 /* ---------------------------------------------------------------------- */
100
101 #ifndef NDEBUG
102 UMF_dump_start ( ) ;
103 init_count = UMF_malloc_count ;
104 DEBUGm4 (("\nUMFPACK numeric: U transpose version\n")) ;
105 #endif
106
107 /* If front_alloc_init negative then allocate that size of front in
108 * UMF_start_front. If alloc_init negative, then allocate that initial
109 * size of Numeric->Memory. */
110
111 relpt = GET_CONTROL (UMFPACK_PIVOT_TOLERANCE,
112 UMFPACK_DEFAULT_PIVOT_TOLERANCE) ;
113 relpt2 = GET_CONTROL (UMFPACK_SYM_PIVOT_TOLERANCE,
114 UMFPACK_DEFAULT_SYM_PIVOT_TOLERANCE) ;
115 alloc_init = GET_CONTROL (UMFPACK_ALLOC_INIT, UMFPACK_DEFAULT_ALLOC_INIT) ;
116 front_alloc_init = GET_CONTROL (UMFPACK_FRONT_ALLOC_INIT,
117 UMFPACK_DEFAULT_FRONT_ALLOC_INIT) ;
118 scale = GET_CONTROL (UMFPACK_SCALE, UMFPACK_DEFAULT_SCALE) ;
119 droptol = GET_CONTROL (UMFPACK_DROPTOL, UMFPACK_DEFAULT_DROPTOL) ;
120
121 relpt = MAX (0.0, MIN (relpt, 1.0)) ;
122 relpt2 = MAX (0.0, MIN (relpt2, 1.0)) ;
123 droptol = MAX (0.0, droptol) ;
124 front_alloc_init = MIN (1.0, front_alloc_init) ;
125
126 if (scale != UMFPACK_SCALE_NONE && scale != UMFPACK_SCALE_MAX)
127 {
128 scale = UMFPACK_DEFAULT_SCALE ;
129 }
130
131 if (User_Info != (double *) NULL)
132 {
133 /* return Info in user's array */
134 Info = User_Info ;
135 /* clear the parts of Info that are set by UMFPACK_numeric */
136 for (i = UMFPACK_NUMERIC_SIZE ; i <= UMFPACK_MAX_FRONT_NCOLS ; i++)
137 {
138 Info [i] = EMPTY ;
139 }
140 for (i = UMFPACK_NUMERIC_DEFRAG ; i < UMFPACK_IR_TAKEN ; i++)
141 {
142 Info [i] = EMPTY ;
143 }
144 }
145 else
146 {
147 /* no Info array passed - use local one instead */
148 Info = Info2 ;
149 for (i = 0 ; i < UMFPACK_INFO ; i++)
150 {
151 Info [i] = EMPTY ;
152 }
153 }
154
155 Symbolic = (SymbolicType *) SymbolicHandle ;
156 Numeric = (NumericType *) NULL ;
157 if (!UMF_valid_symbolic (Symbolic))
158 {
159 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_Symbolic_object ;
160 return (UMFPACK_ERROR_invalid_Symbolic_object) ;
161 }
162
163 /* compute alloc_init automatically for AMD or other symmetric ordering */
164 if (/* Symbolic->ordering == UMFPACK_ORDERING_AMD */ alloc_init >= 0
165 && Symbolic->amd_lunz > 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