1
/* ========================================================================== */
2
/* === UMFPACK_qsymbolic ==================================================== */
3
/* ========================================================================== */
5
/* -------------------------------------------------------------------------- */
6
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A. */
7
/* Davis. All Rights Reserved. See ../README for License. */
8
/* email: davis@cise.ufl.edu CISE Department, Univ. of Florida. */
9
/* web: http://www.cise.ufl.edu/research/sparse/umfpack */
10
/* -------------------------------------------------------------------------- */
13
User-callable. Performs a symbolic factorization.
14
See umfpack_qsymbolic.h and umfpack_symbolic.h for details.
16
Dynamic memory usage: about (3.4nz + 8n + n) integers and n double's as
17
workspace (via UMF_malloc, for a square matrix). All of it is free'd via
18
UMF_free if an error occurs. If successful, the Symbolic object contains
19
12 to 14 objects allocated by UMF_malloc, with a total size of no more
20
than about 13*n integers.
23
#include "umf_internal.h"
24
#include "umf_symbolic_usage.h"
25
#include "umf_colamd.h"
26
#include "umf_set_stats.h"
27
#include "umf_analyze.h"
28
#include "umf_transpose.h"
29
#include "umf_is_permutation.h"
30
#include "umf_malloc.h"
33
#include "umf_singletons.h"
35
typedef struct /* SWType */
37
Int *Front_npivcol ; /* size n_col + 1 */
38
Int *Front_nrows ; /* size n_col */
39
Int *Front_ncols ; /* size n_col */
40
Int *Front_parent ; /* size n_col */
41
Int *Front_cols ; /* size n_col */
42
Int *InFront ; /* size n_row */
43
Int *Ci ; /* size Clen */
44
Int *Cperm1 ; /* size n_col */
45
Int *Rperm1 ; /* size n_row */
46
Int *InvRperm1 ; /* size n_row */
47
Int *Si ; /* size nz */
48
Int *Sp ; /* size n_col + 1 */
49
double *Rs ; /* size n_row */
50
Int *Rperm_2by2 ; /* size n_row */
54
PRIVATE void free_work
61
SymbolicType **Symbolic,
65
/* worst-case usage for SW object */
66
#define SYM_WORK_USAGE(n_col,n_row,Clen) \
67
(DUNITS (Int, Clen) + \
69
4 * DUNITS (Int, n_row) + \
70
4 * DUNITS (Int, n_col) + \
71
2 * DUNITS (Int, n_col + 1) + \
72
DUNITS (double, n_row))
74
/* required size of Ci for code that calls UMF_transpose and UMF_analyze below*/
75
#define UMF_ANALYZE_CLEN(nz,n_row,n_col,nn) \
76
((n_col) + MAX ((nz),(n_col)) + 3*(nn)+1 + (n_col))
78
/* size of an element (in Units), including tuples */
79
#define ELEMENT_SIZE(r,c) \
80
(DGET_ELEMENT_SIZE (r, c) + 1 + (r + c) * UNITS (Tuple, 1))
83
PRIVATE Int init_count ;
86
/* ========================================================================== */
87
/* === do_amd =============================================================== */
88
/* ========================================================================== */
93
const Int Ap [ ], /* size n+1 */
94
const Int Ai [ ], /* size nz = Ap [n] */
95
Int Q [ ], /* output permutation, j = Q [k] */
96
Int Qinv [ ], /* output inverse permutation, Qinv [j] = k */
97
Int Sdeg [ ], /* degree of A+A', from AMD_aat */
98
Int Clen, /* size of Ci */
99
Int Ci [ ], /* size Ci workspace */
100
double amd_Control [ ], /* AMD control parameters */
101
double amd_Info [ ], /* AMD info */
102
SymbolicType *Symbolic, /* Symbolic object */
103
double Info [ ] /* UMFPACK info */
109
Symbolic->amd_dmax = 0 ;
110
Symbolic->amd_lunz = 0 ;
111
Info [UMFPACK_SYMMETRIC_LUNZ] = 0 ;
112
Info [UMFPACK_SYMMETRIC_FLOPS] = 0 ;
113
Info [UMFPACK_SYMMETRIC_DMAX] = 0 ;
114
Info [UMFPACK_SYMMETRIC_NDENSE] = 0 ;
118
AMD_1 (n, Ap, Ai, Q, Qinv, Sdeg, Clen, Ci, amd_Control, amd_Info) ;
120
/* return estimates computed from AMD on PA+PA' */
121
Symbolic->amd_dmax = amd_Info [AMD_DMAX] ;
122
Symbolic->amd_lunz = 2 * amd_Info [AMD_LNZ] + n ;
123
Info [UMFPACK_SYMMETRIC_LUNZ] = Symbolic->amd_lunz ;
124
Info [UMFPACK_SYMMETRIC_FLOPS] = DIV_FLOPS * amd_Info [AMD_NDIV] +
125
MULTSUB_FLOPS * amd_Info [AMD_NMULTSUBS_LU] ;
126
Info [UMFPACK_SYMMETRIC_DMAX] = Symbolic->amd_dmax ;
127
Info [UMFPACK_SYMMETRIC_NDENSE] = amd_Info [AMD_NDENSE] ;
128
Info [UMFPACK_SYMBOLIC_DEFRAG] += amd_Info [AMD_NCMPA] ;
132
/* ========================================================================== */
133
/* === prune_singletons ===================================================== */
134
/* ========================================================================== */
136
/* Create the submatrix after removing the n1 singletons. The matrix has
137
* row and column indices in the range 0 to n_row-n1 and 0 to n_col-n1,
140
PRIVATE Int prune_singletons
160
Int row, k, pp, p, oldcol, newcol, newrow, nzdiag, do_nzdiag ;
163
do_nzdiag = (Ax != (double *) NULL)
165
&& (Az != (double *) NULL)
170
DEBUGm4 (("Prune : S = A (Cperm1 (n1+1:end), Rperm1 (n1+1:end))\n")) ;
171
for (k = 0 ; k < n_row ; k++)
173
ASSERT (Rperm1 [k] >= 0 && Rperm1 [k] < n_row) ;
174
ASSERT (InvRperm1 [Rperm1 [k]] == k) ;
178
/* create the submatrix after removing singletons */
181
for (k = n1 ; k < n_col ; k++)
183
oldcol = Cperm1 [k] ;
185
DEBUG5 (("Prune singletons k "ID" oldcol "ID" newcol "ID": "ID"\n",
186
k, oldcol, newcol, pp)) ;
187
Sp [newcol] = pp ; /* load column pointers */
188
for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
191
DEBUG5 ((" "ID": row "ID, pp, row)) ;
192
ASSERT (row >= 0 && row < n_row) ;
193
newrow = InvRperm1 [row] - n1 ;
194
ASSERT (newrow < n_row - n1) ;
197
DEBUG5 ((" newrow "ID, newrow)) ;
201
/* count the number of truly nonzero entries on the
202
* diagonal of S, excluding entries that are present,
203
* but numerically zero */
204
if (newrow == newcol)
206
/* this is the diagonal entry */
207
if (SCALAR_IS_NONZERO (Ax [p])
209
|| SCALAR_IS_NONZERO (Az [p])
221
Sp [n_col - n1] = pp ;
223
ASSERT (AMD_valid (n_row - n1, n_col - n1, Sp, Si)) ;
227
/* ========================================================================== */
228
/* === combine_ordering ===================================================== */
229
/* ========================================================================== */
231
PRIVATE void combine_ordering
236
Int Cperm_init [ ], /* output permutation */
237
Int Cperm1 [ ], /* singleton and empty column ordering */
238
Int Qinv [ ] /* Qinv from AMD or COLAMD */
241
Int k, oldcol, newcol, knew ;
243
/* combine the singleton ordering with Qinv */
245
for (k = 0 ; k < n_col ; k++)
247
Cperm_init [k] = EMPTY ;
250
for (k = 0 ; k < n1 ; k++)
252
DEBUG1 ((ID" Initial singleton: "ID"\n", k, Cperm1 [k])) ;
253
Cperm_init [k] = Cperm1 [k] ;
255
for (k = n1 ; k < n_col - nempty_col ; k++)
257
/* this is a non-singleton column */
258
oldcol = Cperm1 [k] ; /* user's name for this column */
259
newcol = k - n1 ; /* Qinv's name for this column */
260
knew = Qinv [newcol] ; /* Qinv's ordering for this column */
261
knew += n1 ; /* shift order, after singletons */
262
DEBUG1 ((" k "ID" oldcol "ID" newcol "ID" knew "ID"\n",
263
k, oldcol, newcol, knew)) ;
264
ASSERT (knew >= 0 && knew < n_col - nempty_col) ;
265
ASSERT (Cperm_init [knew] == EMPTY) ;
266
Cperm_init [knew] = oldcol ;
268
for (k = n_col - nempty_col ; k < n_col ; k++)
270
Cperm_init [k] = Cperm1 [k] ;
274
Int *W = (Int *) malloc ((n_col + 1) * sizeof (Int)) ;
275
ASSERT (UMF_is_permutation (Cperm_init, W, n_col, n_col)) ;
282
/* ========================================================================== */
283
/* === UMFPACK_qsymbolic ==================================================== */
284
/* ========================================================================== */
286
GLOBAL Int UMFPACK_qsymbolic
297
void **SymbolicHandle,
298
const double Control [UMFPACK_CONTROL],
299
double User_Info [UMFPACK_INFO]
303
/* ---------------------------------------------------------------------- */
304
/* local variables */
305
/* ---------------------------------------------------------------------- */
307
Int i, nz, j, newj, status, f1, f2, maxnrows, maxncols, nfr, col,
308
nchains, maxrows, maxcols, p, nb, nn, *Chain_start, *Chain_maxrows,
309
*Chain_maxcols, *Front_npivcol, *Ci, Clen, colamd_stats [COLAMD_STATS],
310
fpiv, n_inner, child, parent, *Link, row, *Front_parent,
311
analyze_compactions, k, chain, is_sym, *Si, *Sp, n2, do_UMF_analyze,
312
fpivcol, fallrows, fallcols, *InFront, *F1, snz, *Front_1strow, f1rows,
313
kk, *Cperm_init, *Rperm_init, newrow, *InvRperm1, *Front_leftmostdesc,
314
Clen_analyze, strategy, Clen_amd, fixQ, prefer_diagonal, nzdiag, nzaat,
315
*Wq, *Sdeg, *Fr_npivcol, nempty, *Fr_nrows, *Fr_ncols, *Fr_parent,
316
*Fr_cols, nempty_row, nempty_col, user_auto_strategy, fail, max_rdeg,
317
head_usage, tail_usage, lnz, unz, esize, *Esize, rdeg, *Cdeg, *Rdeg,
318
*Cperm1, *Rperm1, n1, oldcol, newcol, n1c, n1r, *Rperm_2by2, oldrow,
319
dense_row_threshold, tlen, aggressive ;
320
double knobs [COLAMD_KNOBS], flops, f, r, c, *Info, force_fixQ,
321
Info2 [UMFPACK_INFO], drow, dcol, dtail_usage, dlf, duf, dmax_usage,
322
dhead_usage, dlnz, dunz, dmaxfrsize, dClen, dClen_analyze, sym,
323
amd_Info [AMD_INFO], dClen_amd, dr, dc, cr, cc, cp,
324
amd_Control [AMD_CONTROL], stats [2], tol, scale ;
325
SymbolicType *Symbolic ;
326
SWType SWspace, *SW ;
330
init_count = UMF_malloc_count ;
332
"**** Debugging enabled (UMFPACK will be exceedingly slow!) *****************\n"
336
/* ---------------------------------------------------------------------- */
337
/* get the amount of time used by the process so far */
338
/* ---------------------------------------------------------------------- */
340
umfpack_tic (stats) ;
342
/* ---------------------------------------------------------------------- */
343
/* get control settings and check input parameters */
344
/* ---------------------------------------------------------------------- */
346
drow = GET_CONTROL (UMFPACK_DENSE_ROW, UMFPACK_DEFAULT_DENSE_ROW) ;
347
dcol = GET_CONTROL (UMFPACK_DENSE_COL, UMFPACK_DEFAULT_DENSE_COL) ;
348
nb = GET_CONTROL (UMFPACK_BLOCK_SIZE, UMFPACK_DEFAULT_BLOCK_SIZE) ;
349
strategy = GET_CONTROL (UMFPACK_STRATEGY, UMFPACK_DEFAULT_STRATEGY) ;
350
tol = GET_CONTROL (UMFPACK_2BY2_TOLERANCE, UMFPACK_DEFAULT_2BY2_TOLERANCE) ;
351
scale = GET_CONTROL (UMFPACK_SCALE, UMFPACK_DEFAULT_SCALE) ;
352
force_fixQ = GET_CONTROL (UMFPACK_FIXQ, UMFPACK_DEFAULT_FIXQ) ;
353
AMD_defaults (amd_Control) ;
354
amd_Control [AMD_DENSE] =
355
GET_CONTROL (UMFPACK_AMD_DENSE, UMFPACK_DEFAULT_AMD_DENSE) ;
357
(GET_CONTROL (UMFPACK_AGGRESSIVE, UMFPACK_DEFAULT_AGGRESSIVE) != 0) ;
358
amd_Control [AMD_AGGRESSIVE] = aggressive ;
361
nb = MIN (nb, MAXNB) ;
363
if (nb % 2 == 1) nb++ ; /* make sure nb is even */
364
DEBUG0 (("UMFPACK_qsymbolic: nb = "ID" aggressive = "ID"\n", nb,
367
tol = MAX (0.0, MIN (tol, 1.0)) ;
368
if (scale != UMFPACK_SCALE_NONE && scale != UMFPACK_SCALE_MAX)
370
scale = UMFPACK_DEFAULT_SCALE ;
373
if (User_Info != (double *) NULL)
375
/* return Info in user's array */
380
/* no Info array passed - use local one instead */
383
/* clear all of Info */
384
for (i = 0 ; i < UMFPACK_INFO ; i++)
389
nn = MAX (n_row, n_col) ;
390
n_inner = MIN (n_row, n_col) ;
392
Info [UMFPACK_STATUS] = UMFPACK_OK ;
393
Info [UMFPACK_NROW] = n_row ;
394
Info [UMFPACK_NCOL] = n_col ;
395
Info [UMFPACK_SIZE_OF_UNIT] = (double) (sizeof (Unit)) ;
396
Info [UMFPACK_SIZE_OF_INT] = (double) (sizeof (int)) ;
397
Info [UMFPACK_SIZE_OF_LONG] = (double) (sizeof (long)) ;
398
Info [UMFPACK_SIZE_OF_POINTER] = (double) (sizeof (void *)) ;
399
Info [UMFPACK_SIZE_OF_ENTRY] = (double) (sizeof (Entry)) ;
400
Info [UMFPACK_SYMBOLIC_DEFRAG] = 0 ;
402
if (!Ai || !Ap || !SymbolicHandle)
404
Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
405
return (UMFPACK_ERROR_argument_missing) ;
408
*SymbolicHandle = (void *) NULL ;
410
if (n_row <= 0 || n_col <= 0) /* n_row, n_col must be > 0 */
412
Info [UMFPACK_STATUS] = UMFPACK_ERROR_n_nonpositive ;
413
return (UMFPACK_ERROR_n_nonpositive) ;
417
DEBUG0 (("n_row "ID" n_col "ID" nz "ID"\n", n_row, n_col, nz)) ;
418
Info [UMFPACK_NZ] = nz ;
421
Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_matrix ;
422
return (UMFPACK_ERROR_invalid_matrix) ;
425
/* ---------------------------------------------------------------------- */
426
/* get the requested strategy */
427
/* ---------------------------------------------------------------------- */
431
/* if the matrix is rectangular, the only available strategy is
433
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
434
DEBUGm3 (("Rectangular: forcing unsymmetric strategy\n")) ;
437
if (strategy < UMFPACK_STRATEGY_AUTO
438
|| strategy > UMFPACK_STRATEGY_SYMMETRIC)
440
/* unrecognized strategy */
441
strategy = UMFPACK_STRATEGY_AUTO ;
444
user_auto_strategy = (strategy == UMFPACK_STRATEGY_AUTO) ;
446
/* ---------------------------------------------------------------------- */
447
/* determine amount of memory required for UMFPACK_symbolic */
448
/* ---------------------------------------------------------------------- */
450
/* The size of Clen required for UMF_colamd is always larger than */
451
/* UMF_analyze, but the max is included here in case that changes in */
452
/* future versions. */
454
/* This is about 2.2*nz + 9*n_col + 6*n_row, or nz/5 + 13*n_col + 6*n_row,
455
* whichever is bigger. For square matrices, it works out to
456
* 2.2nz + 15n, or nz/5 + 19n, whichever is bigger (typically 2.2nz+15n). */
457
dClen = UMF_COLAMD_RECOMMENDED ((double) nz, (double) n_row,
460
/* This is defined above, as max (nz,n_col) + 3*nn+1 + 2*n_col, where
461
* nn = max (n_row,n_col). It is always smaller than the space required
462
* for colamd or amd. */
463
dClen_analyze = UMF_ANALYZE_CLEN ((double) nz, (double) n_row,
464
(double) n_col, (double) nn) ;
465
dClen = MAX (dClen, dClen_analyze) ;
467
/* The space for AMD can be larger than what's required for colamd: */
468
dClen_amd = 2.4 * (double) nz + 8 * (double) n_inner ;
469
/* additional space for the 2-by-2 strategy */
470
dClen_amd += (double) MAX (nn, nz) ;
471
dClen = MAX (dClen, dClen_amd) ;
473
/* worst case total memory usage for UMFPACK_symbolic (revised below) */
474
Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] =
475
SYM_WORK_USAGE (n_col, n_row, dClen) +
476
UMF_symbolic_usage (n_row, n_col, n_col, n_col, n_col, TRUE) ;
478
if (INT_OVERFLOW (dClen * sizeof (Int)))
480
/* :: int overflow, Clen too large :: */
481
/* Problem is too large for array indexing (Ci [i]) with an Int i. */
482
/* Cannot even analyze the problem to determine upper bounds on */
483
/* memory usage. Need to use the long integer version, umfpack_*l_*. */
484
DEBUGm4 (("out of memory: symbolic int overflow\n")) ;
485
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
486
return (UMFPACK_ERROR_out_of_memory) ;
489
/* repeat the size calculations, in integers */
490
Clen = UMF_COLAMD_RECOMMENDED (nz, n_row, n_col) ;
491
Clen_analyze = UMF_ANALYZE_CLEN (nz, n_row, n_col, nn) ;
492
Clen = MAX (Clen, Clen_analyze) ;
493
Clen_amd = 2.4 * nz + 8 * n_inner ;
494
Clen_amd += MAX (nn, nz) ; /* for Ri, in UMF_2by2 */
495
Clen = MAX (Clen, Clen_amd) ;
497
/* ---------------------------------------------------------------------- */
498
/* allocate the first part of the Symbolic object (header and Cperm_init) */
499
/* ---------------------------------------------------------------------- */
501
/* (1) Five calls to UMF_malloc are made, for a total space of
502
* 2 * (n_row + n_col) + 4 integers + sizeof (SymbolicType).
503
* sizeof (SymbolicType) is a small constant. This space is part of the
504
* Symbolic object and is not freed unless an error occurs. If A is square
505
* then this is about 4*n integers.
508
Symbolic = (SymbolicType *) UMF_malloc (1, sizeof (SymbolicType)) ;
512
/* If we fail here, Symbolic is NULL and thus it won't be */
513
/* dereferenced by UMFPACK_free_symbolic, as called by error ( ). */
514
DEBUGm4 (("out of memory: symbolic object\n")) ;
515
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
516
error (&Symbolic, (SWType *) NULL) ;
517
return (UMFPACK_ERROR_out_of_memory) ;
520
/* We now know that Symbolic has been allocated */
521
Symbolic->valid = 0 ;
522
Symbolic->Chain_start = (Int *) NULL ;
523
Symbolic->Chain_maxrows = (Int *) NULL ;
524
Symbolic->Chain_maxcols = (Int *) NULL ;
525
Symbolic->Front_npivcol = (Int *) NULL ;
526
Symbolic->Front_parent = (Int *) NULL ;
527
Symbolic->Front_1strow = (Int *) NULL ;
528
Symbolic->Front_leftmostdesc = (Int *) NULL ;
529
Symbolic->Esize = (Int *) NULL ;
530
Symbolic->esize = 0 ;
532
Symbolic->Cperm_init = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
533
Symbolic->Rperm_init = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
534
Symbolic->Cdeg = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
535
Symbolic->Rdeg = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
536
Symbolic->Diagonal_map = (Int *) NULL ;
538
Cperm_init = Symbolic->Cperm_init ;
539
Rperm_init = Symbolic->Rperm_init ;
540
Cdeg = Symbolic->Cdeg ;
541
Rdeg = Symbolic->Rdeg ;
543
if (!Cperm_init || !Rperm_init || !Cdeg || !Rdeg)
545
DEBUGm4 (("out of memory: symbolic perm\n")) ;
546
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
547
error (&Symbolic, (SWType *) NULL) ;
548
return (UMFPACK_ERROR_out_of_memory) ;
551
Symbolic->n_row = n_row ;
552
Symbolic->n_col = n_col ;
556
/* ---------------------------------------------------------------------- */
557
/* check user's input permutation */
558
/* ---------------------------------------------------------------------- */
560
if (Quser != (Int *) NULL)
562
/* use Cperm_init as workspace to check input permutation */
563
if (!UMF_is_permutation (Quser, Cperm_init, n_col, n_col))
565
Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_permutation ;
566
error (&Symbolic, (SWType *) NULL) ;
567
return (UMFPACK_ERROR_invalid_permutation) ;
571
/* ---------------------------------------------------------------------- */
572
/* allocate workspace */
573
/* ---------------------------------------------------------------------- */
575
/* (2) Eleven calls to UMF_malloc are made, for workspace of size
576
* Clen + nz + 7*n_col + 2*n_row + 2 integers. Clen is the larger of
577
* MAX (2*nz, 4*n_col) + 8*n_col + 6*n_row + n_col + nz/5 and
578
* 2.4*nz + 8 * MIN (n_row, n_col) + MAX (n_row, n_col, nz)
579
* If A is square and non-singular, then Clen is
580
* MAX (MAX (2*nz, 4*n) + 7*n + nz/5, 3.4*nz) + 8*n
581
* If A has at least 4*n nonzeros then Clen is
582
* MAX (2.2*nz + 7*n, 3.4*nz) + 8*n
583
* If A has at least (7/1.2)*n nonzeros, (about 5.8*n), then Clen is
585
* This space will be free'd when this routine finishes.
587
* Total space thus far is about 3.4nz + 12n integers.
588
* For the double precision, 32-bit integer version, the user's matrix
589
* requires an equivalent space of 3*nz + n integers. So this space is just
590
* slightly larger than the user's input matrix (including the numerical
591
* values themselves).
594
SW = &SWspace ; /* used for UMFPACK_symbolic only */
596
/* Note that SW->Front_* does not include the dummy placeholder front. */
597
/* This space is accounted for by the SYM_WORK_USAGE macro. */
599
/* this is free'd early */
600
SW->Si = (Int *) UMF_malloc (nz, sizeof (Int)) ;
601
SW->Sp = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
602
SW->InvRperm1 = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
603
SW->Cperm1 = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
605
/* this is free'd late */
606
SW->Ci = (Int *) UMF_malloc (Clen, sizeof (Int)) ;
607
SW->Front_npivcol = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
608
SW->Front_nrows = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
609
SW->Front_ncols = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
610
SW->Front_parent = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
611
SW->Front_cols = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
612
SW->Rperm1 = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
613
SW->InFront = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
615
/* this is allocated later, and free'd after Cperm1 but before Ci */
616
SW->Rperm_2by2 = (Int *) NULL ; /* will be nn Int's */
618
/* this is allocated last, and free'd first */
619
SW->Rs = (double *) NULL ; /* will be n_row double's */
622
Fr_npivcol = SW->Front_npivcol ;
623
Fr_nrows = SW->Front_nrows ;
624
Fr_ncols = SW->Front_ncols ;
625
Fr_parent = SW->Front_parent ;
626
Fr_cols = SW->Front_cols ;
627
Cperm1 = SW->Cperm1 ;
628
Rperm1 = SW->Rperm1 ;
631
InvRperm1 = SW->InvRperm1 ;
632
Rperm_2by2 = (Int *) NULL ;
633
InFront = SW->InFront ;
635
if (!Ci || !Fr_npivcol || !Fr_nrows || !Fr_ncols || !Fr_parent || !Fr_cols
636
|| !Cperm1 || !Rperm1 || !Si || !Sp || !InvRperm1 || !InFront)
638
DEBUGm4 (("out of memory: symbolic work\n")) ;
639
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
640
error (&Symbolic, SW) ;
641
return (UMFPACK_ERROR_out_of_memory) ;
644
DEBUG0 (("Symbolic UMF_malloc_count - init_count = "ID"\n",
645
UMF_malloc_count - init_count)) ;
646
ASSERT (UMF_malloc_count == init_count + 17) ;
648
/* ---------------------------------------------------------------------- */
649
/* find the row and column singletons */
650
/* ---------------------------------------------------------------------- */
652
/* [ use first nz + n_row + MAX (n_row, n_col) entries in Ci as workspace,
653
* and use Rperm_init as workspace */
654
ASSERT (Clen >= nz + n_row + MAX (n_row, n_col)) ;
656
status = UMF_singletons (n_row, n_col, Ap, Ai, Quser, Cdeg, Cperm1, Rdeg,
657
Rperm1, InvRperm1, &n1, &n1c, &n1r, &nempty_col, &nempty_row, &is_sym,
658
&max_rdeg, /* workspace: */ Rperm_init, Ci, Ci + nz, Ci + nz + n_row) ;
660
/* ] done using Rperm_init and Ci as workspace */
662
/* InvRperm1 is now the inverse of Rperm1 */
664
if (status != UMFPACK_OK)
666
DEBUGm4 (("matrix invalid: UMF_singletons\n")) ;
667
Info [UMFPACK_STATUS] = status ;
668
error (&Symbolic, SW) ;
671
Info [UMFPACK_NEMPTY_COL] = nempty_col ;
672
Info [UMFPACK_NEMPTY_ROW] = nempty_row ;
673
Info [UMFPACK_NDENSE_COL] = 0 ; /* # dense rows/cols recomputed below */
674
Info [UMFPACK_NDENSE_ROW] = 0 ;
675
Info [UMFPACK_COL_SINGLETONS] = n1c ;
676
Info [UMFPACK_ROW_SINGLETONS] = n1r ;
677
Info [UMFPACK_S_SYMMETRIC] = is_sym ;
679
nempty = MIN (nempty_col, nempty_row) ;
680
Symbolic->nempty_row = nempty_row ;
681
Symbolic->nempty_col = nempty_col ;
683
/* UMF_singletons has verified that the user's input matrix is valid */
684
ASSERT (AMD_valid (n_row, n_col, Ap, Ai)) ;
687
Symbolic->nempty = nempty ;
688
ASSERT (n1 <= n_inner) ;
689
n2 = nn - n1 - nempty ;
691
dense_row_threshold =
692
UMFPACK_DENSE_DEGREE_THRESHOLD (drow, n_col - n1 - nempty_col) ;
693
Symbolic->dense_row_threshold = dense_row_threshold ;
697
/* either the pruned submatrix rectangular, or it is square and
698
* Rperm [n1 .. n-nempty-1] is not the same as Cperm [n1 .. n-nempty-1].
699
* For the auto strategy, switch to the unsymmetric strategy.
700
* Otherwise, if the strategy selected by the user is symmetric or
701
* 2-by-2, then the singletons will be discarded. */
702
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
703
DEBUGm4 (("Strategy: Unsymmetric singletons\n")) ;
706
/* ---------------------------------------------------------------------- */
707
/* determine symmetry, nzdiag, and degrees of S+S' */
708
/* ---------------------------------------------------------------------- */
710
/* S is the matrix obtained after removing singletons
711
* = A (Cperm1 [n1..n_col-nempty_col-1], Rperm1 [n1..n_row-nempty_row-1])
714
Wq = Rperm_init ; /* use Rperm_init as workspace for Wq [ */
715
Sdeg = Cperm_init ; /* use Cperm_init as workspace for Sdeg [ */
719
for (i = 0 ; i < AMD_INFO ; i++)
721
amd_Info [i] = EMPTY ;
724
if (strategy != UMFPACK_STRATEGY_UNSYMMETRIC)
726
/* This also determines the degree of each node in S+S' (Sdeg), which
727
* is needed by the 2-by-2 strategy, the symmetry of S, and the number
728
* of nonzeros on the diagonal of S. */
729
ASSERT (n_row == n_col) ;
730
ASSERT (nempty_row == nempty_col) ;
732
/* get the count of nonzeros on the diagonal of S, excluding explicitly
733
* zero entries. nzdiag = amd_Info [AMD_NZDIAG] counts the zero entries
736
nzdiag = prune_singletons (n1, nn, Ap, Ai, Ax,
740
Cperm1, InvRperm1, Si, Sp
745
nzaat = AMD_aat (n2, Sp, Si, Sdeg, Wq, amd_Info) ;
746
sym = amd_Info [AMD_SYMMETRY] ;
747
Info [UMFPACK_N2] = n2 ;
748
/* nzdiag = amd_Info [AMD_NZDIAG] counts the zero entries of S too */
751
for (k = 0 ; k < n2 ; k++)
753
ASSERT (Sdeg [k] >= 0 && Sdeg [k] < n2) ;
755
ASSERT (Sp [n2] - n2 <= nzaat && nzaat <= 2 * Sp [n2]) ;
756
DEBUG0 (("Explicit zeros: "ID" %g\n", nzdiag, amd_Info [AMD_NZDIAG])) ;
760
/* get statistics from amd_aat, if computed */
761
Symbolic->sym = sym ;
762
Symbolic->nzaat = nzaat ;
763
Symbolic->nzdiag = nzdiag ;
764
Symbolic->amd_dmax = EMPTY ;
766
Info [UMFPACK_PATTERN_SYMMETRY] = sym ;
767
Info [UMFPACK_NZ_A_PLUS_AT] = nzaat ;
768
Info [UMFPACK_NZDIAG] = nzdiag ;
770
/* ---------------------------------------------------------------------- */
771
/* determine the initial strategy based on symmetry and nnz (diag (S)) */
772
/* ---------------------------------------------------------------------- */
774
if (strategy == UMFPACK_STRATEGY_AUTO)
778
/* highly unsymmetric: use the unsymmetric strategy */
779
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
780
DEBUGm4 (("Strategy: select unsymmetric\n")) ;
782
else if (sym >= 0.7 && nzdiag == n2)
784
/* mostly symmetric, zero-free diagonal: use symmetric strategy */
785
strategy = UMFPACK_STRATEGY_SYMMETRIC ;
786
DEBUGm4 (("Strategy: select symmetric\n")) ;
790
/* Evaluate the symmetric 2-by-2 strategy, and select it, or
791
* the unsymmetric strategy if the 2-by-2 strategy doesn't look
793
strategy = UMFPACK_STRATEGY_2BY2 ;
794
DEBUGm4 (("Strategy: try 2-by-2\n")) ;
798
/* ---------------------------------------------------------------------- */
799
/* try the 2-by-2 strategy */
800
/* ---------------------------------------------------------------------- */
802
/* (3) If the 2-by-2 strategy is attempted, additional workspace of size
803
* nn integers and nn double's is allocated, where nn = n_row = n_col.
804
* The real workspace is immediately free'd. The integer workspace of
805
* size nn remains until the end of umfpack_qsymbolic. */
807
/* If the resulting matrix S (Rperm_2by2, :) is too unsymmetric, then the
808
* unsymmetric strategy will be used instead. */
810
if (strategy == UMFPACK_STRATEGY_2BY2)
812
Int *Rp, *Ri, *Blen, *W, nz_papat, nzd2, nweak, unmatched,
816
/* ------------------------------------------------------------------ */
817
/* get workspace for UMF_2by2 */
818
/* ------------------------------------------------------------------ */
820
ASSERT (n_row == n_col && nn == n_row) ;
823
for (k = 0 ; k < n2 ; k++)
825
ASSERT (Sdeg [k] >= 0 && Sdeg [k] < n2) ;
829
/* allocate Rperm_2by2 */
830
SW->Rperm_2by2 = (Int *) UMF_malloc (nn, sizeof (Int)) ;
831
Rperm_2by2 = SW->Rperm_2by2 ;
832
if (Rperm_2by2 == (Int *) NULL)
834
DEBUGm4 (("out of memory: Rperm_2by2\n")) ;
835
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
836
error (&Symbolic, SW) ;
837
return (UMFPACK_ERROR_out_of_memory) ;
840
/* allocate Ri from the tail end of Ci [ */
841
Clen3 = Clen - (MAX (nn, nz) + 1) ;
843
ASSERT (Clen3 >= nz) ; /* space required for UMF_2by2 */
845
/* use Fr_* as workspace for Rp, Blen, and W [ */
850
if (scale != UMFPACK_SCALE_NONE)
852
SW->Rs = (double *) UMF_malloc (nn, sizeof (double)) ;
853
if (SW->Rs == (double *) NULL)
855
DEBUGm4 (("out of memory: scale factors for 2-by-2\n")) ;
856
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
857
error (&Symbolic, SW) ;
858
return (UMFPACK_ERROR_out_of_memory) ;
862
/* ------------------------------------------------------------------ */
863
/* find the 2-by-2 row permutation */
864
/* ------------------------------------------------------------------ */
866
/* find a row permutation Rperm_2by2 such that S (Rperm_2by2, :)
867
* has a healthy diagonal */
869
UMF_2by2 (nn, Ap, Ai, Ax,
877
InvRperm1, n1, nempty, Sdeg, Rperm_2by2, &nweak, &unmatched,
878
Ri, Rp, SW->Rs, Blen, W, Ci, Wq) ;
879
DEBUGm3 (("2by2: nweak "ID" unmatched "ID"\n", nweak, unmatched)) ;
880
Info [UMFPACK_2BY2_NWEAK] = nweak ;
881
Info [UMFPACK_2BY2_UNMATCHED] = unmatched ;
883
SW->Rs = (double *) UMF_free ((void *) SW->Rs) ;
885
/* R = S (Rperm_2by2,:)' */
886
(void) UMF_transpose (n2, n2, Sp, Si, (double *) NULL, Rperm_2by2,
887
(Int *) NULL, 0, Rp, Ri, (double *) NULL, W, FALSE
889
, (double *) NULL, (double *) NULL, FALSE
892
ASSERT (AMD_valid (n2, n2, Rp, Ri)) ;
894
/* contents of Si and Sp no longer needed, but the space is
897
/* ------------------------------------------------------------------ */
898
/* find symmetry of S (Rperm_2by2, :)', and prepare to order with AMD */
899
/* ------------------------------------------------------------------ */
901
for (i = 0 ; i < AMD_INFO ; i++)
903
amd_Info [i] = EMPTY ;
905
nz_papat = AMD_aat (n2, Rp, Ri, Sdeg, Wq, amd_Info) ;
906
sym2 = amd_Info [AMD_SYMMETRY] ;
907
nzd2 = amd_Info [AMD_NZDIAG] ;
909
Info [UMFPACK_2BY2_PATTERN_SYMMETRY] = sym2 ;
910
Info [UMFPACK_2BY2_NZ_PA_PLUS_PAT] = nz_papat ;
911
Info [UMFPACK_2BY2_NZDIAG] = nzd2 ;
913
DEBUG0 (("2by2: sym2 %g nzd2 "ID" n2 "ID"\n", sym2, nzd2, n2)) ;
915
/* ------------------------------------------------------------------ */
916
/* evaluate the 2-by-2 results */
917
/* ------------------------------------------------------------------ */
919
if (user_auto_strategy)
921
if ((sym2 > 1.1 * sym) && (nzd2 > 0.9 * n2))
923
/* 2-by-2 made it much more symmetric */
924
DEBUGm4 (("eval Strategy 2by2: much more symmetric: 2by2\n")) ;
925
strategy = UMFPACK_STRATEGY_2BY2 ;
927
else if (sym2 < 0.7 * sym)
929
/* 2-by-2 made it much more unsymmetric */
930
DEBUGm4 (("eval Strategy 2by2: much more UNsymmetric:unsym\n"));
931
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
933
else if (sym2 < 0.25)
935
DEBUGm4 (("eval Strategy 2by2: is UNsymmetric: unsym\n"));
936
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
938
else if (sym2 >= 0.51)
940
DEBUGm4 (("eval Strategy 2by2: sym2 >= 0.51: 2by2\n")) ;
941
strategy = UMFPACK_STRATEGY_2BY2 ;
943
else if (sym2 >= 0.999 * sym)
945
/* 2-by-2 improved symmetry, or made it only slightly worse */
946
DEBUGm4 (("eval Strategy 2by2: sym2 >= 0.999 sym: 2by2\n")) ;
947
strategy = UMFPACK_STRATEGY_2BY2 ;
951
/* can't decide what to do, so pick the unsymmetric strategy */
952
DEBUGm4 (("eval Strategy 2by2: punt: unsym\n"));
953
strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
957
/* ------------------------------------------------------------------ */
958
/* if the 2-by-2 strategy is selected: */
959
/* ------------------------------------------------------------------ */
961
if (strategy == UMFPACK_STRATEGY_2BY2)
963
if (Quser == (Int *) NULL)
965
/* 2-by-2 strategy is successful */
966
/* compute amd (S) */
967
Int *Qinv = Fr_npivcol ;
968
ASSERT (Clen3 >= (nz_papat + nz_papat/5 + nn) + 7*nn) ;
969
do_amd (n2, Rp, Ri, Wq, Qinv, Sdeg, Clen3, Ci,
970
amd_Control, amd_Info, Symbolic, Info) ;
971
/* combine the singleton ordering and the AMD ordering */
972
combine_ordering (n1, nempty, nn, Cperm_init, Cperm1, Qinv) ;
974
/* fix Rperm_2by2 to reflect A, not S */
975
for (k = 0 ; k < n1 ; k++)
977
oldcol = Cperm1 [k] ;
979
oldrow = Rperm1 [k] ;
980
W [oldcol] = oldrow ;
982
for (k = n1 ; k < nn - nempty ; k++)
984
oldcol = Cperm1 [k] ;
985
i = Rperm_2by2 [k - n1] + n1 ;
986
oldrow = Rperm1 [i] ;
987
W [oldcol] = oldrow ;
989
for (k = nn - nempty ; k < nn ; k++)
991
oldcol = Cperm1 [k] ;
993
oldrow = Rperm1 [k] ;
994
W [oldcol] = oldrow ;
996
for (k = 0 ; k < nn ; k++)
998
Rperm_2by2 [k] = W [k] ;
1001
/* Now, the "diagonal" entry in oldcol (where oldcol is the user's
1002
* name for a column, is the entry in row oldrow (where oldrow is
1003
* the user's name for a row, and oldrow = Rperm_2by2 [oldcol] */
1006
/* Fr_* no longer needed for Rp, Blen, W ] */
1009
/* ---------------------------------------------------------------------- */
1010
/* finalize the strategy, including fixQ and prefer_diagonal */
1011
/* ---------------------------------------------------------------------- */
1013
if (strategy == UMFPACK_STRATEGY_SYMMETRIC)
1015
/* use given Quser or AMD (A+A'), fix Q during factorization,
1016
* prefer diagonal */
1017
DEBUG0 (("\nStrategy: symmetric\n")) ;
1018
ASSERT (n_row == n_col) ;
1019
Symbolic->ordering = UMFPACK_ORDERING_AMD ;
1021
prefer_diagonal = TRUE ;
1023
else if (strategy == UMFPACK_STRATEGY_2BY2)
1025
/* use Q = given Quser or Q = AMD (PA+PA'), fix Q during factorization,
1026
* prefer diagonal, and factorize PAQ, where P is found by UMF_2by2. */
1027
DEBUG0 (("\nStrategy: symmetric 2-by-2\n")) ;
1028
ASSERT (n_row == n_col) ;
1029
Symbolic->ordering = UMFPACK_ORDERING_AMD ;
1031
prefer_diagonal = TRUE ;
1035
/* use given Quser or COLAMD (A), refine Q during factorization,
1036
* no diagonal preference */
1037
ASSERT (strategy == UMFPACK_STRATEGY_UNSYMMETRIC) ;
1038
DEBUG0 (("\nStrategy: unsymmetric\n")) ;
1039
Symbolic->ordering = UMFPACK_ORDERING_COLAMD ;
1041
prefer_diagonal = FALSE ;
1044
if (Quser != (Int *) NULL)
1046
Symbolic->ordering = UMFPACK_ORDERING_GIVEN ;
1052
DEBUG0 (("Force fixQ true\n")) ;
1054
else if (force_fixQ < 0)
1057
DEBUG0 (("Force fixQ false\n")) ;
1060
DEBUG0 (("Strategy: ordering: "ID"\n", Symbolic->ordering)) ;
1061
DEBUG0 (("Strategy: fixQ: "ID"\n", fixQ)) ;
1062
DEBUG0 (("Strategy: prefer diag "ID"\n", prefer_diagonal)) ;
1064
/* get statistics from amd_aat, if computed */
1065
Symbolic->strategy = strategy ;
1066
Symbolic->fixQ = fixQ ;
1067
Symbolic->prefer_diagonal = prefer_diagonal ;
1069
Info [UMFPACK_STRATEGY_USED] = strategy ;
1070
Info [UMFPACK_ORDERING_USED] = Symbolic->ordering ;
1071
Info [UMFPACK_QFIXED] = fixQ ;
1072
Info [UMFPACK_DIAG_PREFERRED] = prefer_diagonal ;
1074
/* ---------------------------------------------------------------------- */
1075
/* get the AMD ordering for the symmetric strategy */
1076
/* ---------------------------------------------------------------------- */
1078
if (strategy == UMFPACK_STRATEGY_SYMMETRIC && Quser == (Int *) NULL)
1080
/* symmetric strategy for a matrix with mostly symmetric pattern */
1081
Int *Qinv = Fr_npivcol ;
1082
ASSERT (n_row == n_col && nn == n_row) ;
1083
ASSERT (Clen >= (nzaat + nzaat/5 + nn) + 7*nn) ;
1084
do_amd (n2, Sp, Si, Wq, Qinv, Sdeg, Clen, Ci,
1085
amd_Control, amd_Info, Symbolic, Info) ;
1086
/* combine the singleton ordering and the AMD ordering */
1087
combine_ordering (n1, nempty, nn, Cperm_init, Cperm1, Qinv) ;
1089
/* Sdeg no longer needed ] */
1090
/* done using Rperm_init as workspace for Wq ] */
1092
/* Contents of Si and Sp no longer needed, but the space is still needed */
1094
/* ---------------------------------------------------------------------- */
1095
/* use the user's input column ordering (already in Cperm1) */
1096
/* ---------------------------------------------------------------------- */
1098
if (Quser != (Int *) NULL)
1100
for (k = 0 ; k < n_col ; k++)
1102
Cperm_init [k] = Cperm1 [k] ;
1106
/* ---------------------------------------------------------------------- */
1107
/* use COLAMD to order the matrix */
1108
/* ---------------------------------------------------------------------- */
1110
if (strategy == UMFPACK_STRATEGY_UNSYMMETRIC && Quser == (Int *) NULL)
1113
/* ------------------------------------------------------------------ */
1114
/* copy the matrix into colamd workspace (colamd destroys its input) */
1115
/* ------------------------------------------------------------------ */
1117
/* C = A (Cperm1 (n1+1:end), Rperm1 (n1+1:end)), where Ci is used as
1118
* the row indices and Cperm_init (on input) is used as the column
1121
(void) prune_singletons (n1, n_col, Ap, Ai,
1126
Cperm1, InvRperm1, Ci, Cperm_init
1132
/* ------------------------------------------------------------------ */
1133
/* set UMF_colamd defaults */
1134
/* ------------------------------------------------------------------ */
1136
UMF_colamd_set_defaults (knobs) ;
1137
knobs [COLAMD_DENSE_ROW] = drow ;
1138
knobs [COLAMD_DENSE_COL] = dcol ;
1139
knobs [COLAMD_AGGRESSIVE] = aggressive ;
1141
/* ------------------------------------------------------------------ */
1142
/* check input matrix and find the initial column pre-ordering */
1143
/* ------------------------------------------------------------------ */
1145
/* NOTE: umf_colamd is not given any original empty rows or columns.
1146
* Those have already been removed via prune_singletons, above. The
1147
* umf_colamd routine has been modified to assume that all rows and
1148
* columns have at least one entry in them. It will break if it is
1149
* given empty rows or columns (an assertion is triggered when running
1153
n_row - n1 - nempty_row,
1154
n_col - n1 - nempty_col,
1155
Clen, Ci, Cperm_init, knobs, colamd_stats,
1156
Fr_npivcol, Fr_nrows, Fr_ncols, Fr_parent, Fr_cols, &nfr,
1158
ASSERT (colamd_stats [COLAMD_EMPTY_ROW] == 0) ;
1159
ASSERT (colamd_stats [COLAMD_EMPTY_COL] == 0) ;
1161
/* # of dense rows will be recomputed below */
1162
Info [UMFPACK_NDENSE_ROW] = colamd_stats [COLAMD_DENSE_ROW] ;
1163
Info [UMFPACK_NDENSE_COL] = colamd_stats [COLAMD_DENSE_COL] ;
1164
Info [UMFPACK_SYMBOLIC_DEFRAG] = colamd_stats [COLAMD_DEFRAG_COUNT] ;
1166
/* re-analyze if any "dense" rows or cols ignored by UMF_colamd */
1168
colamd_stats [COLAMD_DENSE_ROW] > 0 ||
1169
colamd_stats [COLAMD_DENSE_COL] > 0 ;
1171
/* Combine the singleton and colamd ordering into Cperm_init */
1172
/* Note that colamd returns its inverse permutation in Ci */
1173
combine_ordering (n1, nempty_col, n_col, Cperm_init, Cperm1, Ci) ;
1175
/* contents of Ci no longer needed */
1178
for (col = 0 ; col < n_col ; col++)
1180
DEBUG1 (("Cperm_init ["ID"] = "ID"\n", col, Cperm_init[col]));
1182
/* make sure colamd returned a valid permutation */
1183
ASSERT (Cperm_init != (Int *) NULL) ;
1184
ASSERT (UMF_is_permutation (Cperm_init, Ci, n_col, n_col)) ;
1191
/* ------------------------------------------------------------------ */
1192
/* do not call colamd - use input Quser or AMD instead */
1193
/* ------------------------------------------------------------------ */
1195
/* The ordering (Quser or Qamd) is already in Cperm_init */
1196
do_UMF_analyze = TRUE ;
1200
Cperm_init [n_col] = EMPTY ; /* unused in Cperm_init */
1202
/* ---------------------------------------------------------------------- */
1203
/* AMD ordering, if it exists, has been copied into Cperm_init */
1204
/* ---------------------------------------------------------------------- */
1207
DEBUG3 (("Cperm_init column permutation:\n")) ;
1208
ASSERT (UMF_is_permutation (Cperm_init, Ci, n_col, n_col)) ;
1209
for (k = 0 ; k < n_col ; k++)
1211
DEBUG3 ((ID"\n", Cperm_init [k])) ;
1213
/* ensure that empty columns have been placed last in A (:,Cperm_init) */
1214
for (newj = 0 ; newj < n_col ; newj++)
1216
/* empty columns will be last in A (:, Cperm_init (1:n_col)) */
1217
j = Cperm_init [newj] ;
1218
ASSERT (IMPLIES (newj >= n_col-nempty_col, Cdeg [j] == 0)) ;
1219
ASSERT (IMPLIES (newj < n_col-nempty_col, Cdeg [j] > 0)) ;
1223
/* ---------------------------------------------------------------------- */
1224
/* symbolic factorization (unless colamd has already done it) */
1225
/* ---------------------------------------------------------------------- */
1230
Int *W, *Bp, *Bi, *Cperm2, ok, *P, Clen2, bsize, Clen0 ;
1232
/* ------------------------------------------------------------------ */
1233
/* construct column pre-ordered, pruned submatrix */
1234
/* ------------------------------------------------------------------ */
1236
/* S = column form submatrix after removing singletons and applying
1237
* initial column ordering (includes singleton ordering) */
1238
(void) prune_singletons (n1, n_col, Ap, Ai,
1243
Cperm_init, InvRperm1, Si, Sp
1249
/* ------------------------------------------------------------------ */
1250
/* Ci [0 .. Clen-1] holds the following work arrays:
1252
first Clen0 entries empty space, where Clen0 =
1253
Clen - (nn+1 + 2*nn + n_col)
1254
and Clen0 >= nz + n_col
1255
next nn+1 entries Bp [0..nn]
1256
next nn entries Link [0..nn-1]
1257
next nn entries W [0..nn-1]
1258
last n_col entries Cperm2 [0..n_col-1]
1260
We have Clen >= n_col + MAX (nz,n_col) + 3*nn+1 + n_col,
1261
So Clen0 >= 2*n_col as required for AMD_postorder
1262
and Clen0 >= n_col + nz as required
1265
Clen0 = Clen - (nn+1 + 2*nn + n_col) ;
1267
Link = Bp + (nn+1) ;
1270
ASSERT (Cperm2 + n_col == Ci + Clen) ;
1271
ASSERT (Clen0 >= nz + n_col) ;
1272
ASSERT (Clen0 >= 2*n_col) ;
1274
/* ------------------------------------------------------------------ */
1275
/* P = order that rows will be used in UMF_analyze */
1276
/* ------------------------------------------------------------------ */
1278
/* use W to mark rows, and use Link for row permutation P [ [ */
1279
for (row = 0 ; row < n_row - n1 ; row++)
1287
for (col = 0 ; col < n_col - n1 ; col++)
1289
/* empty columns are last in S */
1290
for (p = Sp [col] ; p < Sp [col+1] ; p++)
1295
/* this row has just been seen for the first time */
1302
/* If the matrix has truly empty rows, then P will not be */
1303
/* complete, and visa versa. The matrix is structurally singular. */
1304
nempty_row = n_row - n1 - k ;
1307
/* complete P by putting empty rows last in their natural order, */
1308
/* rather than declaring an error (the matrix is singular) */
1309
for (row = 0 ; row < n_row - n1 ; row++)
1313
/* W [row] = TRUE ; (not required) */
1319
/* contents of W no longer needed ] */
1322
DEBUG3 (("Induced row permutation:\n")) ;
1323
ASSERT (k == n_row - n1) ;
1324
ASSERT (UMF_is_permutation (P, W, n_row - n1, n_row - n1)) ;
1325
for (k = 0 ; k < n_row - n1 ; k++)
1327
DEBUG3 ((ID"\n", P [k])) ;
1331
/* ------------------------------------------------------------------ */
1332
/* B = row-form of the pattern of S (excluding empty columns) */
1333
/* ------------------------------------------------------------------ */
1335
/* Ci [0 .. Clen-1] holds the following work arrays:
1337
first Clen2 entries empty space, must be at least >= n_col
1338
next max (nz,1) Bi [0..max (nz,1)-1]
1339
next nn+1 entries Bp [0..nn]
1340
next nn entries Link [0..nn-1]
1341
next nn entries W [0..nn-1]
1342
last n_col entries Cperm2 [0..n_col-1]
1344
This memory usage is accounted for by the UMF_ANALYZE_CLEN
1349
snz = Sp [n_col - n1] ;
1350
bsize = MAX (snz, 1) ;
1353
ASSERT (Clen2 >= n_col) ;
1355
(void) UMF_transpose (n_row - n1, n_col - n1 - nempty_col,
1356
Sp, Si, (double *) NULL,
1357
P, (Int *) NULL, 0, Bp, Bi, (double *) NULL, W, FALSE
1359
, (double *) NULL, (double *) NULL, FALSE
1363
/* contents of Si and Sp no longer needed */
1365
/* contents of P (same as Link) and W not needed */
1366
/* still need Link and W as work arrays, though ] */
1368
ASSERT (Bp [0] == 0) ;
1369
ASSERT (Bp [n_row - n1] == snz) ;
1371
/* increment Bp to point into Ci, not Bi */
1372
for (i = 0 ; i <= n_row - n1 ; i++)
1376
ASSERT (Bp [0] == Clen0 - bsize) ;
1377
ASSERT (Bp [n_row - n1] <= Clen0) ;
1379
/* Ci [0 .. Clen-1] holds the following work arrays:
1381
first Clen0 entries Ci [0 .. Clen0-1], where the col indices
1382
of B are at the tail end of this part,
1383
and Bp [0] = Clen2 >= n_col. Note
1384
that Clen0 = Clen2 + max (snz,1).
1385
next nn+1 entries Bp [0..nn]
1386
next nn entries Link [0..nn-1]
1387
next nn entries W [0..nn-1]
1388
last n_col entries Cperm2 [0..n_col-1]
1391
/* ------------------------------------------------------------------ */
1393
/* ------------------------------------------------------------------ */
1395
/* only analyze the non-empty, non-singleton part of the matrix */
1397
n_row - n1 - nempty_row,
1398
n_col - n1 - nempty_col,
1399
Ci, Bp, Cperm2, fixQ, W, Link,
1400
Fr_ncols, Fr_nrows, Fr_npivcol,
1401
Fr_parent, &nfr, &analyze_compactions) ;
1404
/* :: internal error in umf_analyze :: */
1405
Info [UMFPACK_STATUS] = UMFPACK_ERROR_internal_error ;
1406
error (&Symbolic, SW) ;
1407
return (UMFPACK_ERROR_internal_error) ;
1409
Info [UMFPACK_SYMBOLIC_DEFRAG] += analyze_compactions ;
1411
/* ------------------------------------------------------------------ */
1412
/* combine the input permutation and UMF_analyze's permutation */
1413
/* ------------------------------------------------------------------ */
1417
/* Cperm2 is the column etree post-ordering */
1418
ASSERT (UMF_is_permutation (Cperm2, W,
1419
n_col-n1-nempty_col, n_col-n1-nempty_col)) ;
1421
/* Note that the empty columns remain at the end of Cperm_init */
1422
for (k = 0 ; k < n_col - n1 - nempty_col ; k++)
1424
W [k] = Cperm_init [n1 + Cperm2 [k]] ;
1427
for (k = 0 ; k < n_col - n1 - nempty_col ; k++)
1429
Cperm_init [n1 + k] = W [k] ;
1433
ASSERT (UMF_is_permutation (Cperm_init, W, n_col, n_col)) ;
1437
/* ---------------------------------------------------------------------- */
1438
/* free some of the workspace */
1439
/* ---------------------------------------------------------------------- */
1441
/* (4) The real workspace, Rs, of size n_row doubles has already been
1442
* free'd. An additional workspace of size nz + n_col+1 + n_col integers
1443
* is now free'd as well. */
1445
SW->Si = (Int *) UMF_free ((void *) SW->Si) ;
1446
SW->Sp = (Int *) UMF_free ((void *) SW->Sp) ;
1447
SW->Cperm1 = (Int *) UMF_free ((void *) SW->Cperm1) ;
1448
ASSERT (SW->Rs == (double *) NULL) ;
1450
/* ---------------------------------------------------------------------- */
1451
/* determine the size of the Symbolic object */
1452
/* ---------------------------------------------------------------------- */
1454
/* ---------------------------------------------------------------------- */
1455
/* determine the size of the Symbolic object */
1456
/* ---------------------------------------------------------------------- */
1459
for (i = 0 ; i < nfr ; i++)
1461
if (Fr_parent [i] != i+1)
1467
Symbolic->nchains = nchains ;
1468
Symbolic->nfr = nfr ;
1470
= (max_rdeg > dense_row_threshold) ? (n_col - n1 - nempty_col) : 0 ;
1472
/* true size of Symbolic object */
1473
Info [UMFPACK_SYMBOLIC_SIZE] = UMF_symbolic_usage (n_row, n_col, nchains,
1474
nfr, Symbolic->esize, prefer_diagonal) ;
1476
/* actual peak memory usage for UMFPACK_symbolic (actual nfr, nchains) */
1477
Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] =
1478
SYM_WORK_USAGE (n_col, n_row, Clen) + Info [UMFPACK_SYMBOLIC_SIZE] ;
1479
Symbolic->peak_sym_usage = Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] ;
1481
DEBUG0 (("Number of fronts: "ID"\n", nfr)) ;
1483
/* ---------------------------------------------------------------------- */
1484
/* allocate the second part of the Symbolic object (Front_*, Chain_*) */
1485
/* ---------------------------------------------------------------------- */
1487
/* (5) UMF_malloc is called 7 or 8 times, for a total space of
1488
* (4*(nfr+1) + 3*(nchains+1) + esize) integers, where nfr is the total
1489
* number of frontal matrices and nchains is the total number of frontal
1490
* matrix chains, and where nchains <= nfr <= n_col. esize is zero if there
1491
* are no dense rows, or n_col-n1-nempty_col otherwise (n1 is the number of
1492
* singletons and nempty_col is the number of empty columns). This space is
1493
* part of the Symbolic object and is not free'd unless an error occurs.
1494
* This is between 7 and about 8n integers when A is square.
1497
/* Note that Symbolic->Front_* does include the dummy placeholder front */
1498
Symbolic->Front_npivcol = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1499
Symbolic->Front_parent = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1500
Symbolic->Front_1strow = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1501
Symbolic->Front_leftmostdesc = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1502
Symbolic->Chain_start = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1503
Symbolic->Chain_maxrows = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1504
Symbolic->Chain_maxcols = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1506
fail = (!Symbolic->Front_npivcol || !Symbolic->Front_parent ||
1507
!Symbolic->Front_1strow || !Symbolic->Front_leftmostdesc ||
1508
!Symbolic->Chain_start || !Symbolic->Chain_maxrows ||
1509
!Symbolic->Chain_maxcols) ;
1511
if (Symbolic->esize > 0)
1513
Symbolic->Esize = (Int *) UMF_malloc (Symbolic->esize, sizeof (Int)) ;
1514
fail = fail || !Symbolic->Esize ;
1519
DEBUGm4 (("out of memory: rest of symbolic object\n")) ;
1520
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
1521
error (&Symbolic, SW) ;
1522
return (UMFPACK_ERROR_out_of_memory) ;
1524
DEBUG0 (("Symbolic UMF_malloc_count - init_count = "ID"\n",
1525
UMF_malloc_count - init_count)) ;
1526
ASSERT (UMF_malloc_count == init_count + 21
1527
+ (SW->Rperm_2by2 != (Int *) NULL)
1528
+ (Symbolic->Esize != (Int *) NULL)) ;
1530
Front_npivcol = Symbolic->Front_npivcol ;
1531
Front_parent = Symbolic->Front_parent ;
1532
Front_1strow = Symbolic->Front_1strow ;
1533
Front_leftmostdesc = Symbolic->Front_leftmostdesc ;
1535
Chain_start = Symbolic->Chain_start ;
1536
Chain_maxrows = Symbolic->Chain_maxrows ;
1537
Chain_maxcols = Symbolic->Chain_maxcols ;
1539
Esize = Symbolic->Esize ;
1541
/* ---------------------------------------------------------------------- */
1542
/* assign rows to fronts */
1543
/* ---------------------------------------------------------------------- */
1545
/* find InFront, unless colamd has already computed it */
1549
DEBUGm4 ((">>>>>>>>>Computing Front_1strow from scratch\n")) ;
1550
/* empty rows go to dummy front nfr */
1551
for (row = 0 ; row < n_row ; row++)
1553
InFront [row] = nfr ;
1555
/* assign the singleton pivot rows to the "empty" front */
1556
for (k = 0 ; k < n1 ; k++)
1559
InFront [row] = EMPTY ;
1561
DEBUG1 (("Front (EMPTY), singleton nrows "ID" ncols "ID"\n", k, k)) ;
1563
for (i = 0 ; i < nfr ; i++)
1565
fpivcol = Fr_npivcol [i] ;
1567
/* for all pivot columns in front i */
1568
for (kk = 0 ; kk < fpivcol ; kk++, newj++)
1570
j = Cperm_init [newj] ;
1571
ASSERT (IMPLIES (newj >= n_col-nempty_col,
1572
Ap [j+1] - Ap [j] == 0));
1573
for (p = Ap [j] ; p < Ap [j+1] ; p++)
1576
if (InFront [row] == nfr)
1578
/* this row belongs to front i */
1579
DEBUG1 ((" Row "ID" in Front "ID"\n", row, i)) ;
1585
Front_1strow [i] = f1rows ;
1586
DEBUG1 ((" Front "ID" has 1strows: "ID" pivcols "ID"\n",
1587
i, f1rows, fpivcol)) ;
1594
/* COLAMD has already computed InFront, but it is not yet
1595
* InFront [row] = front i, where row is an original row. It is
1596
* InFront [k-n1] = i for k in the range n1 to n_row-nempty_row,
1597
* and where row = Rperm1 [k]. Need to permute InFront. Also compute
1598
* # of original rows assembled into each front.
1599
* [ use Ci as workspace */
1600
DEBUGm4 ((">>>>>>>>>Computing Front_1strow from colamd's InFront\n")) ;
1601
for (i = 0 ; i <= nfr ; i++)
1603
Front_1strow [i] = 0 ;
1605
/* assign the singleton pivot rows to "empty" front */
1606
for (k = 0 ; k < n1 ; k++)
1611
/* assign the non-empty rows to the front that assembled them */
1612
for ( ; k < n_row - nempty_row ; k++)
1615
i = InFront [k - n1] ;
1616
ASSERT (i >= EMPTY && i < nfr) ;
1619
Front_1strow [i]++ ;
1621
/* use Ci as permuted version of InFront */
1624
/* empty rows go to the "dummy" front */
1625
for ( ; k < n_row ; k++)
1630
/* permute InFront so that InFront [row] = i if the original row is
1632
for (row = 0 ; row < n_row ; row++)
1634
InFront [row] = Ci [row] ;
1636
/* ] no longer need Ci as workspace */
1640
for (row = 0 ; row < n_row ; row++)
1642
if (InFront [row] == nfr)
1644
DEBUG1 ((" Row "ID" in Dummy Front "ID"\n", row, nfr)) ;
1646
else if (InFront [row] == EMPTY)
1648
DEBUG1 ((" singleton Row "ID"\n", row)) ;
1652
DEBUG1 ((" Row "ID" in Front "ID"\n", row, nfr)) ;
1655
for (i = 0 ; i <= nfr ; i++)
1657
DEBUG1 (("Front "ID" has 1strows: "ID" pivcols "ID"\n",
1658
i, f1rows, fpivcol)) ;
1662
/* ---------------------------------------------------------------------- */
1663
/* copy front information into Symbolic object */
1664
/* ---------------------------------------------------------------------- */
1667
for (i = 0 ; i < nfr ; i++)
1669
fpivcol = Fr_npivcol [i] ;
1670
DEBUG1 (("Front "ID" k "ID" npivcol "ID" nrows "ID" ncols "ID"\n",
1671
i, k, fpivcol, Fr_nrows [i], Fr_ncols [i])) ;
1673
/* copy Front info into Symbolic object from SW */
1674
Front_npivcol [i] = fpivcol ;
1675
Front_parent [i] = Fr_parent [i] ;
1678
/* assign empty columns to dummy placehold front nfr */
1679
DEBUG1 (("Dummy Cols in Front "ID" : "ID"\n", nfr, n_col-k)) ;
1680
Front_npivcol [nfr] = n_col - k ;
1681
Front_parent [nfr] = EMPTY ;
1683
/* ---------------------------------------------------------------------- */
1684
/* find initial row permutation */
1685
/* ---------------------------------------------------------------------- */
1687
/* order the singleton pivot rows */
1688
for (k = 0 ; k < n1 ; k++)
1690
Rperm_init [k] = Rperm1 [k] ;
1693
/* determine the first row in each front (in the new row ordering) */
1694
for (i = 0 ; i < nfr ; i++)
1696
f1rows = Front_1strow [i] ;
1697
DEBUG1 (("Front "ID" : npivcol "ID" parent "ID,
1698
i, Front_npivcol [i], Front_parent [i])) ;
1699
DEBUG1 ((" 1st rows in Front "ID" : "ID"\n", i, f1rows)) ;
1700
Front_1strow [i] = k ;
1704
/* assign empty rows to dummy placehold front nfr */
1705
DEBUG1 (("Rows in Front "ID" (dummy): "ID"\n", nfr, n_row-k)) ;
1706
Front_1strow [nfr] = k ;
1707
DEBUG1 (("nfr "ID" 1strow[nfr] "ID" nrow "ID"\n", nfr, k, n_row)) ;
1709
/* Use Ci as temporary workspace for F1 */
1710
F1 = Ci ; /* [ of size nfr+1 */
1711
ASSERT (Clen >= 2*n_row + nfr+1) ;
1713
for (i = 0 ; i <= nfr ; i++)
1715
F1 [i] = Front_1strow [i] ;
1718
for (row = 0 ; row < n_row ; row++)
1724
ASSERT (newrow >= n1) ;
1725
Rperm_init [newrow] = row ;
1728
Rperm_init [n_row] = EMPTY ; /* unused */
1731
for (k = 0 ; k < n_row ; k++)
1733
DEBUG2 (("Rperm_init ["ID"] = "ID"\n", k, Rperm_init [k])) ;
1737
/* ] done using F1 */
1739
/* ---------------------------------------------------------------------- */
1740
/* find the diagonal map */
1741
/* ---------------------------------------------------------------------- */
1743
/* Rperm_init [newrow] = row gives the row permutation that is implied
1744
* by the column permutation, where "row" is a row index of the original
1745
* matrix A. It is not dependent on the Rperm_2by2 permutation, which
1746
* only redefines the "diagonal". Both are used to construct the
1747
* Diagonal_map. Diagonal_map only needs to be defined for
1748
* k = n1 to nn - nempty, but go ahead and define it for all of
1751
if (prefer_diagonal)
1754
ASSERT (n_row == n_col && nn == n_row) ;
1755
ASSERT (nempty_row == nempty_col && nempty == nempty_row) ;
1757
/* allocate the Diagonal_map */
1758
Symbolic->Diagonal_map = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
1759
Diagonal_map = Symbolic->Diagonal_map ;
1760
if (Diagonal_map == (Int *) NULL)
1762
/* :: out of memory (diagonal map) :: */
1763
DEBUGm4 (("out of memory: Diagonal map\n")) ;
1764
Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
1765
error (&Symbolic, SW) ;
1766
return (UMFPACK_ERROR_out_of_memory) ;
1769
/* use Ci as workspace to compute the inverse of Rperm_init [ */
1770
for (newrow = 0 ; newrow < nn ; newrow++)
1772
oldrow = Rperm_init [newrow] ;
1773
ASSERT (oldrow >= 0 && oldrow < nn) ;
1774
Ci [oldrow] = newrow ;
1776
if (strategy == UMFPACK_STRATEGY_2BY2)
1778
ASSERT (Rperm_2by2 != (Int *) NULL) ;
1779
for (newcol = 0 ; newcol < nn ; newcol++)
1781
oldcol = Cperm_init [newcol] ;
1782
/* 2-by-2 pivoting done in S */
1783
oldrow = Rperm_2by2 [oldcol] ;
1784
newrow = Ci [oldrow] ;
1785
Diagonal_map [newcol] = newrow ;
1790
for (newcol = 0 ; newcol < nn ; newcol++)
1792
oldcol = Cperm_init [newcol] ;
1793
/* no 2-by-2 pivoting in S */
1795
newrow = Ci [oldrow] ;
1796
Diagonal_map [newcol] = newrow ;
1801
DEBUG1 (("\nDiagonal map:\n")) ;
1802
for (newcol = 0 ; newcol < nn ; newcol++)
1804
oldcol = Cperm_init [newcol] ;
1805
DEBUG3 (("oldcol "ID" newcol "ID":\n", oldcol, newcol)) ;
1806
for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
1810
newrow = Ci [oldrow] ;
1811
if (Ax != (double *) NULL
1813
&& Az != (double *) NULL
1817
ASSIGN (aij, Ax [p], Az [p]) ;
1819
if (oldrow == oldcol)
1821
DEBUG2 ((" old diagonal : oldcol "ID" oldrow "ID" ",
1826
if (newrow == Diagonal_map [newcol])
1828
DEBUG2 ((" MAP diagonal : newcol "ID" MAProw "ID" ",
1829
newcol, Diagonal_map [newrow])) ;
1836
/* done using Ci as workspace ] */
1840
/* ---------------------------------------------------------------------- */
1841
/* find the leftmost descendant of each front */
1842
/* ---------------------------------------------------------------------- */
1844
for (i = 0 ; i <= nfr ; i++)
1846
Front_leftmostdesc [i] = EMPTY ;
1849
for (i = 0 ; i < nfr ; i++)
1851
/* start at i and walk up the tree */
1852
DEBUG2 (("Walk up front tree from "ID"\n", i)) ;
1854
while (j != EMPTY && Front_leftmostdesc [j] == EMPTY)
1856
DEBUG3 ((" Leftmost desc of "ID" is "ID"\n", j, i)) ;
1857
Front_leftmostdesc [j] = i ;
1858
j = Front_parent [j] ;
1859
DEBUG3 ((" go to j = "ID"\n", j)) ;
1863
/* ---------------------------------------------------------------------- */
1864
/* find the frontal matrix chains and max frontal matrix sizes */
1865
/* ---------------------------------------------------------------------- */
1867
maxnrows = 1 ; /* max # rows in any front */
1868
maxncols = 1 ; /* max # cols in any front */
1869
dmaxfrsize = 1 ; /* max frontal matrix size */
1871
/* start the first chain */
1872
nchains = 0 ; /* number of chains */
1873
Chain_start [0] = 0 ; /* front 0 starts a new chain */
1874
maxrows = 1 ; /* max # rows for any front in current chain */
1875
maxcols = 1 ; /* max # cols for any front in current chain */
1876
DEBUG1 (("Constructing chains:\n")) ;
1878
for (i = 0 ; i < nfr ; i++)
1880
/* get frontal matrix info */
1881
fpivcol = Front_npivcol [i] ; /* # candidate pivot columns */
1882
fallrows = Fr_nrows [i] ; /* all rows (not just Schur comp) */
1883
fallcols = Fr_ncols [i] ; /* all cols (not just Schur comp) */
1884
parent = Front_parent [i] ; /* parent in column etree */
1885
fpiv = MIN (fpivcol, fallrows) ; /* # pivot rows and cols */
1886
maxrows = MAX (maxrows, fallrows) ;
1887
maxcols = MAX (maxcols, fallcols) ;
1889
DEBUG1 (("Front: "ID", pivcol "ID", "ID"-by-"ID" parent "ID
1890
", npiv "ID" Chain: maxrows "ID" maxcols "ID"\n", i, fpivcol,
1891
fallrows, fallcols, parent, fpiv, maxrows, maxcols)) ;
1895
/* this is the end of a chain */
1897
DEBUG1 (("\nEnd of chain "ID"\n", nchains)) ;
1899
/* make sure maxrows is an odd number */
1900
ASSERT (maxrows >= 0) ;
1901
if (maxrows % 2 == 0) maxrows++ ;
1903
DEBUG1 (("Chain maxrows "ID" maxcols "ID"\n", maxrows, maxcols)) ;
1905
Chain_maxrows [nchains] = maxrows ;
1906
Chain_maxcols [nchains] = maxcols ;
1908
/* keep track of the maximum front size for all chains */
1910
/* for Info only: */
1911
s = (double) maxrows * (double) maxcols ;
1912
dmaxfrsize = MAX (dmaxfrsize, s) ;
1914
/* for the subsequent numerical factorization */
1915
maxnrows = MAX (maxnrows, maxrows) ;
1916
maxncols = MAX (maxncols, maxcols) ;
1918
DEBUG1 (("Chain dmaxfrsize %g\n\n", dmaxfrsize)) ;
1920
/* start the next chain */
1922
Chain_start [nchains] = i+1 ;
1928
/* for Info only: */
1929
dmaxfrsize = ceil (dmaxfrsize) ;
1930
DEBUGm1 (("dmaxfrsize %30.20g Int_MAX "ID"\n", dmaxfrsize, Int_MAX)) ;
1931
ASSERT (Symbolic->nchains == nchains) ;
1933
/* For allocating objects in umfpack_numeric (does not include all possible
1934
* pivots, particularly pivots from prior fronts in the chain. Need to add
1935
* nb to these to get the # of columns in the L block, for example. This
1936
* is the largest row dimension and largest column dimension of any frontal
1937
* matrix. maxnrows is always odd. */
1938
Symbolic->maxnrows = maxnrows ;
1939
Symbolic->maxncols = maxncols ;
1940
DEBUGm3 (("maxnrows "ID" maxncols "ID"\n", maxnrows, maxncols)) ;
1942
/* ---------------------------------------------------------------------- */
1943
/* find the initial element sizes */
1944
/* ---------------------------------------------------------------------- */
1946
if (max_rdeg > dense_row_threshold)
1948
/* there are one or more dense rows in the input matrix */
1949
/* count the number of dense rows in each column */
1950
/* use Ci as workspace for inverse of Rperm_init [ */
1951
ASSERT (Esize != (Int *) NULL) ;
1952
for (newrow = 0 ; newrow < n_row ; newrow++)
1954
oldrow = Rperm_init [newrow] ;
1955
ASSERT (oldrow >= 0 && oldrow < nn) ;
1956
Ci [oldrow] = newrow ;
1958
for (col = n1 ; col < n_col - nempty_col ; col++)
1960
oldcol = Cperm_init [col] ;
1961
esize = Cdeg [oldcol] ;
1962
ASSERT (esize > 0) ;
1963
for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
1966
newrow = Ci [oldrow] ;
1967
if (newrow >= n1 && Rdeg [oldrow] > dense_row_threshold)
1972
ASSERT (esize >= 0) ;
1973
Esize [col - n1] = esize ;
1975
/* done using Ci as workspace ] */
1978
/* If there are no dense rows, then Esize [col-n1] is identical to
1979
* Cdeg [col], once Cdeg is permuted below */
1981
/* ---------------------------------------------------------------------- */
1982
/* permute Cdeg and Rdeg according to initial column and row permutation */
1983
/* ---------------------------------------------------------------------- */
1985
/* use Ci as workspace [ */
1986
for (k = 0 ; k < n_col ; k++)
1988
Ci [k] = Cdeg [Cperm_init [k]] ;
1990
for (k = 0 ; k < n_col ; k++)
1994
for (k = 0 ; k < n_row ; k++)
1996
Ci [k] = Rdeg [Rperm_init [k]] ;
1998
for (k = 0 ; k < n_row ; k++)
2002
/* done using Ci as workspace ] */
2004
/* ---------------------------------------------------------------------- */
2005
/* simulate UMF_kernel_init */
2006
/* ---------------------------------------------------------------------- */
2008
/* count elements and tuples at tail, LU factors of singletons, and
2009
* head and tail markers */
2011
dlnz = n_inner ; /* upper limit of nz in L (incl diag) */
2012
dunz = dlnz ; /* upper limit of nz in U (incl diag) */
2022
/* allocate the Rpi and Rpx workspace for UMF_kernel_init (incl. headers) */
2023
tail_usage += UNITS (Int *, n_row+1) + UNITS (Entry *, n_row+1) + 2 ;
2024
dtail_usage += DUNITS (Int *, n_row+1) + DUNITS (Entry *, n_row+1) + 2 ;
2025
DEBUG1 (("Symbolic usage after Rpi/Rpx allocation: head "ID" tail "ID"\n",
2026
head_usage, tail_usage)) ;
2028
/* LU factors for singletons, at the head of memory */
2029
for (k = 0 ; k < n1 ; k++)
2031
lnz = Cdeg [k] - 1 ;
2032
unz = Rdeg [k] - 1 ;
2035
DEBUG1 (("singleton k "ID" pivrow "ID" pivcol "ID" lnz "ID" unz "ID"\n",
2036
k, Rperm_init [k], Cperm_init [k], lnz, unz)) ;
2037
head_usage += UNITS (Int, lnz) + UNITS (Entry, lnz)
2038
+ UNITS (Int, unz) + UNITS (Entry, unz) ;
2039
dhead_usage += DUNITS (Int, lnz) + DUNITS (Entry, lnz)
2040
+ DUNITS (Int, unz) + DUNITS (Entry, unz) ;
2042
DEBUG1 (("Symbolic init head usage: "ID" for LU singletons\n",head_usage)) ;
2044
/* column elements: */
2045
for (k = n1 ; k < n_col - nempty_col; k++)
2047
esize = Esize ? Esize [k-n1] : Cdeg [k] ;
2048
DEBUG2 ((" esize: "ID"\n", esize)) ;
2049
ASSERT (esize >= 0) ;
2052
tail_usage += GET_ELEMENT_SIZE (esize, 1) + 1 ;
2053
dtail_usage += DGET_ELEMENT_SIZE (esize, 1) + 1 ;
2057
/* dense row elements */
2060
Int nrow_elements = 0 ;
2061
for (k = n1 ; k < n_row - nempty_row ; k++)
2064
if (rdeg > dense_row_threshold)
2066
tail_usage += GET_ELEMENT_SIZE (1, rdeg) + 1 ;
2067
dtail_usage += GET_ELEMENT_SIZE (1, rdeg) + 1 ;
2071
Info [UMFPACK_NDENSE_ROW] = nrow_elements ;
2074
DEBUG1 (("Symbolic usage: "ID" = head "ID" + tail "ID" after els\n",
2075
head_usage + tail_usage, head_usage, tail_usage)) ;
2077
/* compute the tuple lengths */
2081
for (row = n1 ; row < n_row ; row++)
2084
tlen = (rdeg > dense_row_threshold) ? 1 : rdeg ;
2085
tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2086
dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2089
for (col = n1 ; col < n_col - nempty_col ; col++)
2091
/* tlen is 1 plus the number of dense rows in this column */
2092
esize = Esize [col - n1] ;
2093
tlen = (esize > 0) + (Cdeg [col] - esize) ;
2094
tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2095
dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2097
for ( ; col < n_col ; col++)
2099
tail_usage += 1 + UNITS (Tuple, TUPLES (0)) ;
2100
dtail_usage += 1 + DUNITS (Tuple, TUPLES (0)) ;
2106
for (row = n1 ; row < n_row ; row++)
2109
tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2110
dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2113
for (col = n1 ; col < n_col ; col++)
2115
tail_usage += 1 + UNITS (Tuple, TUPLES (1)) ;
2116
dtail_usage += 1 + DUNITS (Tuple, TUPLES (1)) ;
2120
Symbolic->num_mem_init_usage = head_usage + tail_usage ;
2121
DEBUG1 (("Symbolic usage: "ID" = head "ID" + tail "ID" final\n",
2122
Symbolic->num_mem_init_usage, head_usage, tail_usage)) ;
2124
ASSERT (UMF_is_permutation (Rperm_init, Ci, n_row, n_row)) ;
2126
/* initial head and tail usage in Numeric->Memory */
2127
dmax_usage = dhead_usage + dtail_usage ;
2128
dmax_usage = MAX (Symbolic->num_mem_init_usage, ceil (dmax_usage)) ;
2129
Info [UMFPACK_VARIABLE_INIT_ESTIMATE] = dmax_usage ;
2131
/* In case Symbolic->num_mem_init_usage overflows, keep as a double, too */
2132
Symbolic->dnum_mem_init_usage = dmax_usage ;
2134
/* free the Rpi and Rpx workspace */
2135
tail_usage -= UNITS (Int *, n_row+1) + UNITS (Entry *, n_row+1) ;
2136
dtail_usage -= DUNITS (Int *, n_row+1) + DUNITS (Entry *, n_row+1) ;
2138
/* ---------------------------------------------------------------------- */
2139
/* simulate UMF_kernel, assuming unsymmetric pivoting */
2140
/* ---------------------------------------------------------------------- */
2142
/* Use Ci as temporary workspace for link lists [ */
2144
for (i = 0 ; i < nfr ; i++)
2149
flops = 0 ; /* flop count upper bound */
2151
for (chain = 0 ; chain < nchains ; chain++)
2154
f1 = Chain_start [chain] ;
2155
f2 = Chain_start [chain+1] - 1 ;
2157
/* allocate frontal matrix working array (C, L, and U) */
2158
dr = Chain_maxrows [chain] ;
2159
dc = Chain_maxcols [chain] ;
2161
nb*nb /* LU is nb-by-nb */
2162
+ dr*nb /* L is dr-by-nb */
2163
+ nb*dc /* U is nb-by-dc, stored by rows */
2164
+ dr*dc ; /* C is dr by dc */
2165
dtail_usage += DUNITS (Entry, fsize) ;
2166
dmax_usage = MAX (dmax_usage, dhead_usage + dtail_usage) ;
2168
for (i = f1 ; i <= f2 ; i++)
2171
/* get frontal matrix info */
2172
fpivcol = Front_npivcol [i] ; /* # candidate pivot columns */
2173
fallrows = Fr_nrows [i] ; /* all rows (not just Schur comp*/
2174
fallcols = Fr_ncols [i] ; /* all cols (not just Schur comp*/
2175
parent = Front_parent [i] ; /* parent in column etree */
2176
fpiv = MIN (fpivcol, fallrows) ; /* # pivot rows and cols */
2178
r = fallrows - fpiv ; /* # rows in Schur comp. */
2179
c = fallcols - fpiv ; /* # cols in Schur comp. */
2181
/* assemble all children of front i in column etree */
2182
for (child = Link [i] ; child != EMPTY ; child = Link [child])
2184
ASSERT (child >= 0 && child < i) ;
2185
ASSERT (Front_parent [child] == i) ;
2186
/* free the child element and remove it from tuple lists */
2187
cp = MIN (Front_npivcol [child], Fr_nrows [child]) ;
2188
cr = Fr_nrows [child] - cp ;
2189
cc = Fr_ncols [child] - cp ;
2190
ASSERT (cp >= 0 && cr >= 0 && cc >= 0) ;
2191
dtail_usage -= ELEMENT_SIZE (cr, cc) ;
2195
/* The flop count computed here is "canonical". */
2197
/* factorize the frontal matrix */
2198
flops += DIV_FLOPS * (f*r + (f-1)*f/2) /* scale pivot columns */
2199
/* f outer products: */
2200
+ MULTSUB_FLOPS * (f*r*c + (r+c)*(f-1)*f/2 + (f-1)*f*(2*f-1)/6);
2202
/* count nonzeros and memory usage in double precision */
2203
dlf = (f*f-f)/2 + f*r ; /* nz in L below diagonal */
2204
duf = (f*f-f)/2 + f*c ; /* nz in U above diagonal */
2208
/* store f columns of L and f rows of U */
2210
DUNITS (Entry, dlf + duf) /* numerical values (excl diag) */
2211
+ DUNITS (Int, r + c + f) ; /* indices (compressed) */
2213
if (parent != EMPTY)
2215
/* create new element and place in tuple lists */
2216
dtail_usage += ELEMENT_SIZE (r, c) ;
2218
/* place in link list of parent */
2219
Link [i] = Link [parent] ;
2223
/* keep track of peak Numeric->Memory usage */
2224
dmax_usage = MAX (dmax_usage, dhead_usage + dtail_usage) ;
2228
/* free the current frontal matrix */
2229
dtail_usage -= DUNITS (Entry, fsize) ;
2232
dhead_usage = ceil (dhead_usage) ;
2233
dmax_usage = ceil (dmax_usage) ;
2234
Symbolic->num_mem_size_est = dhead_usage ;
2235
Symbolic->num_mem_usage_est = dmax_usage ;
2236
Symbolic->lunz_bound = dlnz + dunz - n_inner ;
2238
/* ] done using Ci as workspace for Link array */
2240
/* ---------------------------------------------------------------------- */
2241
/* estimate total memory usage in UMFPACK_numeric */
2242
/* ---------------------------------------------------------------------- */
2247
dmax_usage, /* estimated peak size of Numeric->Memory */
2248
dhead_usage, /* estimated final size of Numeric->Memory */
2249
flops, /* estimated "true flops" */
2250
dlnz, /* estimated nz in L */
2251
dunz, /* estimated nz in U */
2252
dmaxfrsize, /* estimated largest front size */
2253
(double) n_col, /* worst case Numeric->Upattern size */
2254
(double) n_inner, /* max possible pivots to be found */
2255
(double) maxnrows, /* estimated largest #rows in front */
2256
(double) maxncols, /* estimated largest #cols in front */
2257
TRUE, /* assume scaling is to be performed */
2261
/* ---------------------------------------------------------------------- */
2264
for (i = 0 ; i < nchains ; i++)
2266
DEBUG2 (("Chain "ID" start "ID" end "ID" maxrows "ID" maxcols "ID"\n",
2267
i, Chain_start [i], Chain_start [i+1] - 1,
2268
Chain_maxrows [i], Chain_maxcols [i])) ;
2269
UMF_dump_chain (Chain_start [i], Fr_parent, Fr_npivcol, Fr_nrows,
2273
for (i = 0 ; i < nfr ; i++)
2275
fpivcol = MAX (fpivcol, Front_npivcol [i]) ;
2277
DEBUG0 (("Max pivot cols in any front: "ID"\n", fpivcol)) ;
2278
DEBUG1 (("Largest front: maxnrows "ID" maxncols "ID" dmaxfrsize %g\n",
2279
maxnrows, maxncols, dmaxfrsize)) ;
2282
/* ---------------------------------------------------------------------- */
2283
/* UMFPACK_symbolic was successful, return the object handle */
2284
/* ---------------------------------------------------------------------- */
2286
Symbolic->valid = SYMBOLIC_VALID ;
2287
*SymbolicHandle = (void *) Symbolic ;
2289
/* ---------------------------------------------------------------------- */
2290
/* free workspace */
2291
/* ---------------------------------------------------------------------- */
2293
/* (6) The last of the workspace is free'd. The final Symbolic object
2294
* consists of 12 to 14 allocated objects. Its final total size is lies
2295
* roughly between 4*n and 13*n for a square matrix, which is all that is
2296
* left of the memory allocated by this routine. If an error occurs, the
2297
* entire Symbolic object is free'd when this routine returns (the error
2298
* return routine, below).
2303
DEBUG0 (("(3)Symbolic UMF_malloc_count - init_count = "ID"\n",
2304
UMF_malloc_count - init_count)) ;
2305
ASSERT (UMF_malloc_count == init_count + 12
2306
+ (Symbolic->Esize != (Int *) NULL)
2307
+ (Symbolic->Diagonal_map != (Int *) NULL)) ;
2309
/* ---------------------------------------------------------------------- */
2310
/* get the time used by UMFPACK_*symbolic */
2311
/* ---------------------------------------------------------------------- */
2313
umfpack_toc (stats) ;
2314
Info [UMFPACK_SYMBOLIC_WALLTIME] = stats [0] ;
2315
Info [UMFPACK_SYMBOLIC_TIME] = stats [1] ;
2317
return (UMFPACK_OK) ;
2321
/* ========================================================================== */
2322
/* === free_work ============================================================ */
2323
/* ========================================================================== */
2325
PRIVATE void free_work
2332
SW->Rperm_2by2 = (Int *) UMF_free ((void *) SW->Rperm_2by2) ;
2333
SW->InvRperm1 = (Int *) UMF_free ((void *) SW->InvRperm1) ;
2334
SW->Rs = (double *) UMF_free ((void *) SW->Rs) ;
2335
SW->Si = (Int *) UMF_free ((void *) SW->Si) ;
2336
SW->Sp = (Int *) UMF_free ((void *) SW->Sp) ;
2337
SW->Ci = (Int *) UMF_free ((void *) SW->Ci) ;
2338
SW->Front_npivcol = (Int *) UMF_free ((void *) SW->Front_npivcol);
2339
SW->Front_nrows = (Int *) UMF_free ((void *) SW->Front_nrows) ;
2340
SW->Front_ncols = (Int *) UMF_free ((void *) SW->Front_ncols) ;
2341
SW->Front_parent = (Int *) UMF_free ((void *) SW->Front_parent) ;
2342
SW->Front_cols = (Int *) UMF_free ((void *) SW->Front_cols) ;
2343
SW->Cperm1 = (Int *) UMF_free ((void *) SW->Cperm1) ;
2344
SW->Rperm1 = (Int *) UMF_free ((void *) SW->Rperm1) ;
2345
SW->InFront = (Int *) UMF_free ((void *) SW->InFront) ;
2350
/* ========================================================================== */
2351
/* === error ================================================================ */
2352
/* ========================================================================== */
2354
/* Error return from UMFPACK_symbolic. Free all allocated memory. */
2358
SymbolicType **Symbolic,
2364
UMFPACK_free_symbolic ((void **) Symbolic) ;
2365
ASSERT (UMF_malloc_count == init_count) ;
added
removed
Lib/sandbox/pysparse/umfpack/umfpack_qsymbolic.c
