1
/* ========================================================================== */
2
/* === UMFPACK_numeric ====================================================== */
3
/* ========================================================================== */
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
/* -------------------------------------------------------------------------- */
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
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.
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"
28
#include "umf_realloc.h"
31
PRIVATE Int init_count ;
34
PRIVATE Int work_alloc
37
SymbolicType *Symbolic
40
PRIVATE void free_work
45
PRIVATE Int numeric_alloc
47
NumericType **NumericHandle,
48
SymbolicType *Symbolic,
55
NumericType **Numeric,
60
/* ========================================================================== */
61
/* === UMFPACK_numeric ====================================================== */
62
/* ========================================================================== */
64
GLOBAL Int UMFPACK_numeric
74
const double Control [UMFPACK_CONTROL],
75
double User_Info [UMFPACK_INFO]
79
/* ---------------------------------------------------------------------- */
81
/* ---------------------------------------------------------------------- */
83
double Info2 [UMFPACK_INFO], alloc_init, relpt, relpt2, droptol,
84
front_alloc_init, stats [2] ;
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 ;
92
/* ---------------------------------------------------------------------- */
93
/* get the amount of time used by the process so far */
94
/* ---------------------------------------------------------------------- */
98
/* ---------------------------------------------------------------------- */
99
/* initialize and check inputs */
100
/* ---------------------------------------------------------------------- */
104
init_count = UMF_malloc_count ;
105
DEBUGm4 (("\nUMFPACK numeric: U transpose version\n")) ;
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. */
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) ;
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) ;
127
if (scale != UMFPACK_SCALE_NONE && scale != UMFPACK_SCALE_MAX)
129
scale = UMFPACK_DEFAULT_SCALE ;
132
if (User_Info != (double *) NULL)
134
/* return Info in user's array */
136
/* clear the parts of Info that are set by UMFPACK_numeric */
137
for (i = UMFPACK_NUMERIC_SIZE ; i <= UMFPACK_MAX_FRONT_NCOLS ; i++)
141
for (i = UMFPACK_NUMERIC_DEFRAG ; i < UMFPACK_IR_TAKEN ; i++)
148
/* no Info array passed - use local one instead */
150
for (i = 0 ; i < UMFPACK_INFO ; i++)
156
Symbolic = (SymbolicType *) SymbolicHandle ;
157
Numeric = (NumericType *) NULL ;
158
if (!UMF_valid_symbolic (Symbolic))
160
Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_Symbolic_object ;
161
return (UMFPACK_ERROR_invalid_Symbolic_object) ;
164
/* compute alloc_init automatically for AMD ordering */
165
if (Symbolic->ordering == UMFPACK_ORDERING_AMD && alloc_init >= 0)
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 ;
172
n_row = Symbolic->n_row ;
173
n_col = Symbolic->n_col ;
174
n_inner = MIN (n_row, n_col) ;
176
/* check for integer overflow in Numeric->Memory minimum size */
177
if (INT_OVERFLOW (Symbolic->dnum_mem_init_usage * sizeof (Unit)))
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) ;
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)) ;
192
if (!Ap || !Ai || !Ax || !NumericHandle)
194
Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
195
return (UMFPACK_ERROR_argument_missing) ;
198
Info [UMFPACK_NZ] = Ap [n_col] ;
199
*NumericHandle = (void *) NULL ;
201
/* ---------------------------------------------------------------------- */
202
/* allocate the Work object */
203
/* ---------------------------------------------------------------------- */
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.
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 ;
223
if (!work_alloc (Work, Symbolic))
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) ;
230
ASSERT (UMF_malloc_count == init_count + 16 + 2*Symbolic->prefer_diagonal) ;
232
/* ---------------------------------------------------------------------- */
233
/* allocate Numeric object */
234
/* ---------------------------------------------------------------------- */
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.
247
if (!numeric_alloc (&Numeric, Symbolic, alloc_init, scale))
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) ;
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))) ;
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 ;
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)) ;
272
/* ---------------------------------------------------------------------- */
273
/* scale and factorize */
274
/* ---------------------------------------------------------------------- */
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.
288
DEBUG0 (("Calling umf_kernel\n")) ;
289
status = UMF_kernel (Ap, Ai, Ax,
293
Numeric, Work, Symbolic) ;
295
Info [UMFPACK_STATUS] = status ;
296
if (status < UMFPACK_OK)
298
/* out of memory, or pattern has changed */
299
error (&Numeric, Work) ;
303
Info [UMFPACK_FORCED_UPDATES] = Work->nforced ;
304
Info [UMFPACK_VARIABLE_INIT] = Numeric->init_usage ;
305
if (Symbolic->prefer_diagonal)
307
Info [UMFPACK_NOFF_DIAG] = Work->noff_diagonal ;
310
DEBUG0 (("malloc: init_count "ID" UMF_malloc_count "ID"\n",
311
init_count, UMF_malloc_count)) ;
313
npiv = Numeric->npiv ; /* = n_inner for nonsingular matrices */
314
ulen = Numeric->ulen ; /* = 0 for square nonsingular matrices */
316
/* ---------------------------------------------------------------------- */
317
/* free Work object */
318
/* ---------------------------------------------------------------------- */
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.
327
DEBUG0 (("malloc: init_count "ID" UMF_malloc_count "ID"\n",
328
init_count, UMF_malloc_count)) ;
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))) ;
338
/* ---------------------------------------------------------------------- */
339
/* reduce Lpos, Lilen, Lip, Upos, Uilen and Uip to size npiv+1 */
340
/* ---------------------------------------------------------------------- */
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.
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)) ;
354
Numeric->Lpos = inew ;
356
inew = (Int *) UMF_realloc (Numeric->Uilen, npiv+1, sizeof (Int)) ;
359
Numeric->Uilen = inew ;
361
inew = (Int *) UMF_realloc (Numeric->Uip, npiv+1, sizeof (Int)) ;
364
Numeric->Uip = inew ;
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)) ;
374
Numeric->Upos = inew ;
376
inew = (Int *) UMF_realloc (Numeric->Lilen, npiv+1, sizeof (Int)) ;
379
Numeric->Lilen = inew ;
381
inew = (Int *) UMF_realloc (Numeric->Lip, npiv+1, sizeof (Int)) ;
384
Numeric->Lip = inew ;
388
/* ---------------------------------------------------------------------- */
389
/* reduce Numeric->Upattern from size n_col+1 to size ulen+1 */
390
/* ---------------------------------------------------------------------- */
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
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)
400
inew = (Int *) UMF_realloc (Numeric->Upattern, ulen+1, sizeof (Int)) ;
403
Numeric->Upattern = inew ;
407
/* ---------------------------------------------------------------------- */
408
/* reduce Numeric->Memory to hold just the LU factors at the head */
409
/* ---------------------------------------------------------------------- */
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
416
newsize = Numeric->ihead ;
417
if (newsize < Numeric->size)
419
mnew = (Unit *) UMF_realloc (Numeric->Memory, newsize, sizeof (Unit)) ;
422
/* realloc succeeded (how can it fail since the size is reduced?) */
423
Numeric->Memory = mnew ;
424
Numeric->size = newsize ;
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) */
433
/* ---------------------------------------------------------------------- */
434
/* report the results and return the Numeric object */
435
/* ---------------------------------------------------------------------- */
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,
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 +
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 ;
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) ;
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))
479
/* rcond is zero if there is any zero or NaN on the diagonal */
480
Numeric->rcond = 0.0 ;
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 ;
488
Info [UMFPACK_UMIN] = Numeric->min_udiag ;
489
Info [UMFPACK_UMAX] = Numeric->max_udiag ;
490
Info [UMFPACK_RCOND] = Numeric->rcond ;
492
if (Numeric->nnzpiv < n_inner
493
|| SCALAR_IS_ZERO (Numeric->rcond) || SCALAR_IS_NAN (Numeric->rcond))
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 ;
503
Numeric->valid = NUMERIC_VALID ;
504
*NumericHandle = (void *) Numeric ;
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 */
511
/* ---------------------------------------------------------------------- */
512
/* get the time used by UMFPACK_numeric */
513
/* ---------------------------------------------------------------------- */
515
umfpack_toc (stats) ;
516
Info [UMFPACK_NUMERIC_WALLTIME] = stats [0] ;
517
Info [UMFPACK_NUMERIC_TIME] = stats [1] ;
519
/* return UMFPACK_OK or UMFPACK_WARNING_singular_matrix */
525
/* ========================================================================== */
526
/* === numeric_alloc ======================================================== */
527
/* ========================================================================== */
529
/* Allocate the Numeric object */
531
PRIVATE Int numeric_alloc
533
NumericType **NumericHandle,
534
SymbolicType *Symbolic,
539
double nsize, bsize ;
540
Int n_row, n_col, n_inner, min_usage, trying ;
541
NumericType *Numeric ;
543
DEBUG0 (("numeric alloc:\n")) ;
545
n_row = Symbolic->n_row ;
546
n_col = Symbolic->n_col ;
547
n_inner = MIN (n_row, n_col) ;
548
*NumericHandle = (NumericType *) NULL ;
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)) ;
556
return (FALSE) ; /* out of memory */
559
*NumericHandle = Numeric ;
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)) ;
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)
577
DEBUG0 (("Allocating scale factors\n")) ;
578
Numeric->Rs = (double *) UMF_malloc (n_row, sizeof (double)) ;
582
DEBUG0 (("No scale factors allocated (R = I)\n")) ;
583
Numeric->Rs = (double *) NULL ;
586
Numeric->Memory = (Unit *) NULL ;
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 */
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))
599
return (FALSE) ; /* out of memory */
602
/* ---------------------------------------------------------------------- */
603
/* allocate initial Numeric->Memory for LU factors and elements */
604
/* ---------------------------------------------------------------------- */
608
/* -alloc_init is the exact size to initially allocate */
609
nsize = -alloc_init ;
613
/* alloc_init is a ratio of the upper bound memory usage */
614
nsize = (alloc_init * Symbolic->num_mem_usage_est) + 1 ;
616
min_usage = Symbolic->num_mem_init_usage ;
618
/* Numeric->Memory must be large enough for UMF_kernel_init */
619
nsize = MAX (min_usage, nsize) ;
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) ;
627
Numeric->size = (Int) nsize ;
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)) ;
632
/* allocates 1 object: */
633
/* keep trying until successful, or memory request is too small */
637
Numeric->Memory = (Unit *) UMF_malloc (Numeric->size, sizeof (Unit)) ;
640
DEBUG0 (("Successful Numeric->size: "ID"\n", Numeric->size)) ;
643
/* too much, reduce the request (but not below the minimum) */
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) ;
651
return (FALSE) ; /* we failed to allocate Numeric->Memory */
655
/* ========================================================================== */
656
/* === work_alloc =========================================================== */
657
/* ========================================================================== */
659
/* Allocate the Work object. Return TRUE if successful. */
661
PRIVATE Int work_alloc
664
SymbolicType *Symbolic
667
Int n_row, n_col, nn, maxnrows, maxncols, nfr, ok, maxnrc, n1 ;
669
n_row = Work->n_row ;
670
n_col = Work->n_col ;
671
nn = MAX (n_row, n_col) ;
674
ASSERT (n1 <= n_row && n1 <= n_col) ;
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) ;
682
DEBUG0 (("work alloc: maxnrows+nb "ID" maxncols+nb "ID"\n",
683
maxnrows, maxncols)) ;
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)) ;
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) ;
709
/* 2 allocations: accounted for in UMF_set_stats (work_usage) */
710
if (Symbolic->prefer_diagonal)
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 ;
718
/* no diagonal map needed for rectangular matrices */
719
Work->Diagonal_map = (Int *) NULL ;
720
Work->Diagonal_imap = (Int *) NULL ;
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 ;
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 ;
733
DEBUG0 (("work alloc done.\n")) ;
738
/* ========================================================================== */
739
/* === free_work ============================================================ */
740
/* ========================================================================== */
742
PRIVATE void free_work
747
DEBUG0 (("work free:\n")) ;
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) ;
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) ;
773
DEBUG0 (("work free done.\n")) ;
777
/* ========================================================================== */
778
/* === error ================================================================ */
779
/* ========================================================================== */
781
/* Error return from UMFPACK_numeric. Free all allocated memory. */
785
NumericType **Numeric,
790
UMFPACK_free_numeric ((void **) Numeric) ;
791
ASSERT (UMF_malloc_count == init_count) ;
added
removed
contrib/umfpack/UMFPACK/Source/umfpack_numeric.c
