1 /* ========================================================================== */
2 /* === UMFPACK_qsymbolic ==================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* -------------------------------------------------------------------------- */
10
11 /*
12 User-callable. Performs a symbolic factorization.
13 See umfpack_qsymbolic.h and umfpack_symbolic.h for details.
14
15 Dynamic memory usage: about (3.4nz + 8n + n) integers and n double's as
16 workspace (via UMF_malloc, for a square matrix). All of it is free'd via
17 UMF_free if an error occurs. If successful, the Symbolic object contains
18 12 to 14 objects allocated by UMF_malloc, with a total size of no more
19 than about 13*n integers.
20 */
21
22 #include "umf_internal.h"
23 #include "umf_symbolic_usage.h"
24 #include "umf_colamd.h"
25 #include "umf_set_stats.h"
26 #include "umf_analyze.h"
27 #include "umf_transpose.h"
28 #include "umf_is_permutation.h"
29 #include "umf_malloc.h"
30 #include "umf_free.h"
31 #include "umf_2by2.h"
32 #include "umf_singletons.h"
33
34 typedef struct /* SWType */
35 {
36 Int *Front_npivcol ; /* size n_col + 1 */
37 Int *Front_nrows ; /* size n_col */
38 Int *Front_ncols ; /* size n_col */
39 Int *Front_parent ; /* size n_col */
40 Int *Front_cols ; /* size n_col */
41 Int *InFront ; /* size n_row */
42 Int *Ci ; /* size Clen */
43 Int *Cperm1 ; /* size n_col */
44 Int *Rperm1 ; /* size n_row */
45 Int *InvRperm1 ; /* size n_row */
46 Int *Si ; /* size nz */
47 Int *Sp ; /* size n_col + 1 */
48 double *Rs ; /* size n_row */
49 Int *Rperm_2by2 ; /* size n_row */
50
51 } SWType ;
52
53 PRIVATE void free_work
54 (
55 SWType *SW
56 ) ;
57
58 PRIVATE void error
59 (
60 SymbolicType **Symbolic,
61 SWType *SW
62 ) ;
63
64 /* worst-case usage for SW object */
65 #define SYM_WORK_USAGE(n_col,n_row,Clen) \
66 (DUNITS (Int, Clen) + \
67 DUNITS (Int, nz) + \
68 4 * DUNITS (Int, n_row) + \
69 4 * DUNITS (Int, n_col) + \
70 2 * DUNITS (Int, n_col + 1) + \
71 DUNITS (double, n_row))
72
73 /* required size of Ci for code that calls UMF_transpose and UMF_analyze below*/
74 #define UMF_ANALYZE_CLEN(nz,n_row,n_col,nn) \
75 ((n_col) + MAX ((nz),(n_col)) + 3*(nn)+1 + (n_col))
76
77 /* size of an element (in Units), including tuples */
78 #define ELEMENT_SIZE(r,c) \
79 (DGET_ELEMENT_SIZE (r, c) + 1 + (r + c) * UNITS (Tuple, 1))
80
81 #ifndef NDEBUG
82 PRIVATE Int init_count ;
83 #endif
84
85 /* ========================================================================== */
86 /* === do_amd =============================================================== */
87 /* ========================================================================== */
88
do_amd(Int n,const Int Ap[],const Int Ai[],Int Q[],Int Qinv[],Int Sdeg[],Int Clen,Int Ci[],double amd_Control[],double amd_Info[],SymbolicType * Symbolic,double Info[])89 PRIVATE void do_amd
90 (
91 Int n,
92 const Int Ap [ ], /* size n+1 */
93 const Int Ai [ ], /* size nz = Ap [n] */
94 Int Q [ ], /* output permutation, j = Q [k] */
95 Int Qinv [ ], /* output inverse permutation, Qinv [j] = k */
96 Int Sdeg [ ], /* degree of A+A', from AMD_aat */
97 Int Clen, /* size of Ci */
98 Int Ci [ ], /* size Ci workspace */
99 double amd_Control [ ], /* AMD control parameters */
100 double amd_Info [ ], /* AMD info */
101 SymbolicType *Symbolic, /* Symbolic object */
102 double Info [ ] /* UMFPACK info */
103 )
104 {
105
106 if (n == 0)
107 {
108 Symbolic->amd_dmax = 0 ;
109 Symbolic->amd_lunz = 0 ;
110 Info [UMFPACK_SYMMETRIC_LUNZ] = 0 ;
111 Info [UMFPACK_SYMMETRIC_FLOPS] = 0 ;
112 Info [UMFPACK_SYMMETRIC_DMAX] = 0 ;
113 Info [UMFPACK_SYMMETRIC_NDENSE] = 0 ;
114 }
115 else
116 {
117 AMD_1 (n, Ap, Ai, Q, Qinv, Sdeg, Clen, Ci, amd_Control, amd_Info) ;
118
119 /* return estimates computed from AMD on PA+PA' */
120 Symbolic->amd_dmax = amd_Info [AMD_DMAX] ;
121 Symbolic->amd_lunz = 2 * amd_Info [AMD_LNZ] + n ;
122 Info [UMFPACK_SYMMETRIC_LUNZ] = Symbolic->amd_lunz ;
123 Info [UMFPACK_SYMMETRIC_FLOPS] = DIV_FLOPS * amd_Info [AMD_NDIV] +
124 MULTSUB_FLOPS * amd_Info [AMD_NMULTSUBS_LU] ;
125 Info [UMFPACK_SYMMETRIC_DMAX] = Symbolic->amd_dmax ;
126 Info [UMFPACK_SYMMETRIC_NDENSE] = amd_Info [AMD_NDENSE] ;
127 Info [UMFPACK_SYMBOLIC_DEFRAG] += amd_Info [AMD_NCMPA] ;
128 }
129 }
130
131 /* ========================================================================== */
132 /* === prune_singletons ===================================================== */
133 /* ========================================================================== */
134
135 /* Create the submatrix after removing the n1 singletons. The matrix has
136 * row and column indices in the range 0 to n_row-n1 and 0 to n_col-n1,
137 * respectively. */
138
prune_singletons(Int n1,Int n_col,const Int Ap[],const Int Ai[],const double Ax[],const double Az[],Int Cperm1[],Int InvRperm1[],Int Si[],Int Sp[],Int Rperm1[],Int n_row)139 PRIVATE Int prune_singletons
140 (
141 Int n1,
142 Int n_col,
143 const Int Ap [ ],
144 const Int Ai [ ],
145 const double Ax [ ],
146 #ifdef COMPLEX
147 const double Az [ ],
148 #endif
149 Int Cperm1 [ ],
150 Int InvRperm1 [ ],
151 Int Si [ ],
152 Int Sp [ ]
153 #ifndef NDEBUG
154 , Int Rperm1 [ ]
155 , Int n_row
156 #endif
157 )
158 {
159 Int row, k, pp, p, oldcol, newcol, newrow, nzdiag, do_nzdiag ;
160 #ifdef COMPLEX
161 Int split = SPLIT (Az) ;
162 #endif
163
164 nzdiag = 0 ;
165 do_nzdiag = (Ax != (double *) NULL) ;
166
167 #ifndef NDEBUG
168 DEBUGm4 (("Prune : S = A (Cperm1 (n1+1:end), Rperm1 (n1+1:end))\n")) ;
169 for (k = 0 ; k < n_row ; k++)
170 {
171 ASSERT (Rperm1 [k] >= 0 && Rperm1 [k] < n_row) ;
172 ASSERT (InvRperm1 [Rperm1 [k]] == k) ;
173 }
174 #endif
175
176 /* create the submatrix after removing singletons */
177
178 pp = 0 ;
179 for (k = n1 ; k < n_col ; k++)
180 {
181 oldcol = Cperm1 [k] ;
182 newcol = k - n1 ;
183 DEBUG5 (("Prune singletons k " ID " oldcol " ID " newcol " ID ": " ID "\n",
184 k, oldcol, newcol, pp)) ;
185 Sp [newcol] = pp ; /* load column pointers */
186 for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
187 {
188 row = Ai [p] ;
189 DEBUG5 ((" " ID ": row " ID, pp, row)) ;
190 ASSERT (row >= 0 && row < n_row) ;
191 newrow = InvRperm1 [row] - n1 ;
192 ASSERT (newrow < n_row - n1) ;
193 if (newrow >= 0)
194 {
195 DEBUG5 ((" newrow " ID, newrow)) ;
196 Si [pp++] = newrow ;
197 if (do_nzdiag)
198 {
199 /* count the number of truly nonzero entries on the
200 * diagonal of S, excluding entries that are present,
201 * but numerically zero */
202 if (newrow == newcol)
203 {
204 /* this is the diagonal entry */
205 #ifdef COMPLEX
206 if (split)
207 {
208 if (SCALAR_IS_NONZERO (Ax [p]) ||
209 SCALAR_IS_NONZERO (Az [p]))
210 {
211 nzdiag++ ;
212 }
213 }
214 else
215 {
216 if (SCALAR_IS_NONZERO (Ax [2*p ]) ||
217 SCALAR_IS_NONZERO (Ax [2*p+1]))
218 {
219 nzdiag++ ;
220 }
221 }
222 #else
223 if (SCALAR_IS_NONZERO (Ax [p]))
224 {
225 nzdiag++ ;
226 }
227 #endif
228 }
229 }
230 }
231 DEBUG5 (("\n")) ;
232 }
233 }
234 Sp [n_col - n1] = pp ;
235
236 return (nzdiag) ;
237 }
238
239 /* ========================================================================== */
240 /* === combine_ordering ===================================================== */
241 /* ========================================================================== */
242
combine_ordering(Int n1,Int nempty_col,Int n_col,Int Cperm_init[],Int Cperm1[],Int Qinv[])243 PRIVATE void combine_ordering
244 (
245 Int n1,
246 Int nempty_col,
247 Int n_col,
248 Int Cperm_init [ ], /* output permutation */
249 Int Cperm1 [ ], /* singleton and empty column ordering */
250 Int Qinv [ ] /* Qinv from AMD or COLAMD */
251 )
252 {
253 Int k, oldcol, newcol, knew ;
254
255 /* combine the singleton ordering with Qinv */
256 #ifndef NDEBUG
257 for (k = 0 ; k < n_col ; k++)
258 {
259 Cperm_init [k] = EMPTY ;
260 }
261 #endif
262 for (k = 0 ; k < n1 ; k++)
263 {
264 DEBUG1 ((ID " Initial singleton: " ID "\n", k, Cperm1 [k])) ;
265 Cperm_init [k] = Cperm1 [k] ;
266 }
267 for (k = n1 ; k < n_col - nempty_col ; k++)
268 {
269 /* this is a non-singleton column */
270 oldcol = Cperm1 [k] ; /* user's name for this column */
271 newcol = k - n1 ; /* Qinv's name for this column */
272 knew = Qinv [newcol] ; /* Qinv's ordering for this column */
273 knew += n1 ; /* shift order, after singletons */
274 DEBUG1 ((" k " ID " oldcol " ID " newcol " ID " knew " ID "\n",
275 k, oldcol, newcol, knew)) ;
276 ASSERT (knew >= 0 && knew < n_col - nempty_col) ;
277 ASSERT (Cperm_init [knew] == EMPTY) ;
278 Cperm_init [knew] = oldcol ;
279 }
280 for (k = n_col - nempty_col ; k < n_col ; k++)
281 {
282 Cperm_init [k] = Cperm1 [k] ;
283 }
284 #ifndef NDEBUG
285 {
286 Int *W = (Int *) malloc ((n_col + 1) * sizeof (Int)) ;
287 ASSERT (UMF_is_permutation (Cperm_init, W, n_col, n_col)) ;
288 free (W) ;
289 }
290 #endif
291
292 }
293
294 /* ========================================================================== */
295 /* === UMFPACK_qsymbolic ==================================================== */
296 /* ========================================================================== */
297
UMFPACK_qsymbolic(Int n_row,Int n_col,const Int Ap[],const Int Ai[],const double Ax[],const double Az[],const Int Quser[],void ** SymbolicHandle,const double Control[UMFPACK_CONTROL],double User_Info[UMFPACK_INFO])298 GLOBAL Int UMFPACK_qsymbolic
299 (
300 Int n_row,
301 Int n_col,
302 const Int Ap [ ],
303 const Int Ai [ ],
304 const double Ax [ ],
305 #ifdef COMPLEX
306 const double Az [ ],
307 #endif
308 const Int Quser [ ],
309 void **SymbolicHandle,
310 const double Control [UMFPACK_CONTROL],
311 double User_Info [UMFPACK_INFO]
312 )
313 {
314
315 /* ---------------------------------------------------------------------- */
316 /* local variables */
317 /* ---------------------------------------------------------------------- */
318
319 double knobs [COLAMD_KNOBS], flops, f, r, c, force_fixQ,
320 Info2 [UMFPACK_INFO], drow, dcol, dtail_usage, dlf, duf, dmax_usage,
321 dhead_usage, dlnz, dunz, dmaxfrsize, dClen, dClen_analyze, sym,
322 amd_Info [AMD_INFO], dClen_amd, dr, dc, cr, cc, cp,
323 amd_Control [AMD_CONTROL], stats [2], tol ;
324 double *Info ;
325 Int i, nz, j, newj, status, f1, f2, maxnrows, maxncols, nfr, col,
326 nchains, maxrows, maxcols, p, nb, nn, *Chain_start, *Chain_maxrows,
327 *Chain_maxcols, *Front_npivcol, *Ci, Clen, colamd_stats [COLAMD_STATS],
328 fpiv, n_inner, child, parent, *Link, row, *Front_parent,
329 analyze_compactions, k, chain, is_sym, *Si, *Sp, n2, do_UMF_analyze,
330 fpivcol, fallrows, fallcols, *InFront, *F1, snz, *Front_1strow, f1rows,
331 kk, *Cperm_init, *Rperm_init, newrow, *InvRperm1, *Front_leftmostdesc,
332 Clen_analyze, strategy, Clen_amd, fixQ, prefer_diagonal, nzdiag, nzaat,
333 *Wq, *Sdeg, *Fr_npivcol, nempty, *Fr_nrows, *Fr_ncols, *Fr_parent,
334 *Fr_cols, nempty_row, nempty_col, user_auto_strategy, fail, max_rdeg,
335 head_usage, tail_usage, lnz, unz, esize, *Esize, rdeg, *Cdeg, *Rdeg,
336 *Cperm1, *Rperm1, n1, oldcol, newcol, n1c, n1r, *Rperm_2by2, oldrow,
337 dense_row_threshold, tlen, aggressive, scale, *Rp, *Ri ;
338
339 SymbolicType *Symbolic ;
340 SWType SWspace, *SW ;
341
342 #ifndef NDEBUG
343 UMF_dump_start ( ) ;
344 init_count = UMF_malloc_count ;
345 PRINTF ((
346 "**** Debugging enabled (UMFPACK will be exceedingly slow!) *****************\n"
347 )) ;
348 #endif
349
350 /* ---------------------------------------------------------------------- */
351 /* get the amount of time used by the process so far */
352 /* ---------------------------------------------------------------------- */
353
354 umfpack_tic (stats) ;
355
356 /* ---------------------------------------------------------------------- */
357 /* get control settings and check input parameters */
358 /* ---------------------------------------------------------------------- */
359
360 drow = GET_CONTROL (UMFPACK_DENSE_ROW, UMFPACK_DEFAULT_DENSE_ROW) ;
361 dcol = GET_CONTROL (UMFPACK_DENSE_COL, UMFPACK_DEFAULT_DENSE_COL) ;
362 nb = GET_CONTROL (UMFPACK_BLOCK_SIZE, UMFPACK_DEFAULT_BLOCK_SIZE) ;
363 strategy = GET_CONTROL (UMFPACK_STRATEGY, UMFPACK_DEFAULT_STRATEGY) ;
364 tol = GET_CONTROL (UMFPACK_2BY2_TOLERANCE, UMFPACK_DEFAULT_2BY2_TOLERANCE) ;
365 scale = GET_CONTROL (UMFPACK_SCALE, UMFPACK_DEFAULT_SCALE) ;
366 force_fixQ = GET_CONTROL (UMFPACK_FIXQ, UMFPACK_DEFAULT_FIXQ) ;
367 AMD_defaults (amd_Control) ;
368 amd_Control [AMD_DENSE] =
369 GET_CONTROL (UMFPACK_AMD_DENSE, UMFPACK_DEFAULT_AMD_DENSE) ;
370 aggressive =
371 (GET_CONTROL (UMFPACK_AGGRESSIVE, UMFPACK_DEFAULT_AGGRESSIVE) != 0) ;
372 amd_Control [AMD_AGGRESSIVE] = aggressive ;
373
374 nb = MAX (2, nb) ;
375 nb = MIN (nb, MAXNB) ;
376 ASSERT (nb >= 0) ;
377 if (nb % 2 == 1) nb++ ; /* make sure nb is even */
378 DEBUG0 (("UMFPACK_qsymbolic: nb = " ID " aggressive = " ID "\n", nb,
379 aggressive)) ;
380
381 tol = MAX (0.0, MIN (tol, 1.0)) ;
382 if (scale != UMFPACK_SCALE_NONE && scale != UMFPACK_SCALE_MAX)
383 {
384 scale = UMFPACK_DEFAULT_SCALE ;
385 }
386
387 if (User_Info != (double *) NULL)
388 {
389 /* return Info in user's array */
390 Info = User_Info ;
391 }
392 else
393 {
394 /* no Info array passed - use local one instead */
395 Info = Info2 ;
396 }
397 /* clear all of Info */
398 for (i = 0 ; i < UMFPACK_INFO ; i++)
399 {
400 Info [i] = EMPTY ;
401 }
402
403 nn = MAX (n_row, n_col) ;
404 n_inner = MIN (n_row, n_col) ;
405
406 Info [UMFPACK_STATUS] = UMFPACK_OK ;
407 Info [UMFPACK_NROW] = n_row ;
408 Info [UMFPACK_NCOL] = n_col ;
409 Info [UMFPACK_SIZE_OF_UNIT] = (double) (sizeof (Unit)) ;
410 Info [UMFPACK_SIZE_OF_INT] = (double) (sizeof (int)) ;
411 Info [UMFPACK_SIZE_OF_LONG] = (double) (sizeof (UF_long)) ;
412 Info [UMFPACK_SIZE_OF_POINTER] = (double) (sizeof (void *)) ;
413 Info [UMFPACK_SIZE_OF_ENTRY] = (double) (sizeof (Entry)) ;
414 Info [UMFPACK_SYMBOLIC_DEFRAG] = 0 ;
415
416 if (!Ai || !Ap || !SymbolicHandle)
417 {
418 Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
419 return (UMFPACK_ERROR_argument_missing) ;
420 }
421
422 *SymbolicHandle = (void *) NULL ;
423
424 if (n_row <= 0 || n_col <= 0) /* n_row, n_col must be > 0 */
425 {
426 Info [UMFPACK_STATUS] = UMFPACK_ERROR_n_nonpositive ;
427 return (UMFPACK_ERROR_n_nonpositive) ;
428 }
429
430 nz = Ap [n_col] ;
431 DEBUG0 (("n_row " ID " n_col " ID " nz " ID "\n", n_row, n_col, nz)) ;
432 Info [UMFPACK_NZ] = nz ;
433 if (nz < 0)
434 {
435 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_matrix ;
436 return (UMFPACK_ERROR_invalid_matrix) ;
437 }
438
439 /* ---------------------------------------------------------------------- */
440 /* get the requested strategy */
441 /* ---------------------------------------------------------------------- */
442
443 if (n_row != n_col)
444 {
445 /* if the matrix is rectangular, the only available strategy is
446 * unsymmetric */
447 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
448 DEBUGm3 (("Rectangular: forcing unsymmetric strategy\n")) ;
449 }
450
451 if (strategy < UMFPACK_STRATEGY_AUTO
452 || strategy > UMFPACK_STRATEGY_SYMMETRIC)
453 {
454 /* unrecognized strategy */
455 strategy = UMFPACK_STRATEGY_AUTO ;
456 }
457
458 if (Quser != (Int *) NULL)
459 {
460 /* when the user provides Q, only symmetric and unsymmetric strategies
461 * are available */
462 if (strategy == UMFPACK_STRATEGY_2BY2)
463 {
464 strategy = UMFPACK_STRATEGY_SYMMETRIC ;
465 }
466 if (strategy != UMFPACK_STRATEGY_SYMMETRIC)
467 {
468 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
469 }
470 }
471
472 user_auto_strategy = (strategy == UMFPACK_STRATEGY_AUTO) ;
473
474 /* ---------------------------------------------------------------------- */
475 /* determine amount of memory required for UMFPACK_symbolic */
476 /* ---------------------------------------------------------------------- */
477
478 /* The size of Clen required for UMF_colamd is always larger than */
479 /* UMF_analyze, but the max is included here in case that changes in */
480 /* future versions. */
481
482 /* This is about 2.2*nz + 9*n_col + 6*n_row, or nz/5 + 13*n_col + 6*n_row,
483 * whichever is bigger. For square matrices, it works out to
484 * 2.2nz + 15n, or nz/5 + 19n, whichever is bigger (typically 2.2nz+15n). */
485 dClen = UMF_COLAMD_RECOMMENDED ((double) nz, (double) n_row,
486 (double) n_col) ;
487
488 /* This is defined above, as max (nz,n_col) + 3*nn+1 + 2*n_col, where
489 * nn = max (n_row,n_col). It is always smaller than the space required
490 * for colamd or amd. */
491 dClen_analyze = UMF_ANALYZE_CLEN ((double) nz, (double) n_row,
492 (double) n_col, (double) nn) ;
493 dClen = MAX (dClen, dClen_analyze) ;
494
495 /* The space for AMD can be larger than what's required for colamd: */
496 dClen_amd = 2.4 * (double) nz + 8 * (double) n_inner ;
497 /* additional space for the 2-by-2 strategy */
498 dClen_amd += (double) MAX (nn, nz) ;
499 dClen = MAX (dClen, dClen_amd) ;
500
501 /* worst case total memory usage for UMFPACK_symbolic (revised below) */
502 Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] =
503 SYM_WORK_USAGE (n_col, n_row, dClen) +
504 UMF_symbolic_usage (n_row, n_col, n_col, n_col, n_col, TRUE) ;
505
506 if (INT_OVERFLOW (dClen * sizeof (Int)))
507 {
508 /* :: int overflow, Clen too large :: */
509 /* Problem is too large for array indexing (Ci [i]) with an Int i. */
510 /* Cannot even analyze the problem to determine upper bounds on */
511 /* memory usage. Need to use the UF_long version, umfpack_*l_*. */
512 DEBUGm4 (("out of memory: symbolic int overflow\n")) ;
513 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
514 return (UMFPACK_ERROR_out_of_memory) ;
515 }
516
517 /* repeat the size calculations, in integers */
518 Clen = UMF_COLAMD_RECOMMENDED (nz, n_row, n_col) ;
519 Clen_analyze = UMF_ANALYZE_CLEN (nz, n_row, n_col, nn) ;
520 Clen = MAX (Clen, Clen_analyze) ;
521 Clen_amd = 2.4 * nz + 8 * n_inner ;
522 Clen_amd += MAX (nn, nz) ; /* for Ri, in UMF_2by2 */
523 Clen = MAX (Clen, Clen_amd) ;
524
525 /* ---------------------------------------------------------------------- */
526 /* allocate the first part of the Symbolic object (header and Cperm_init) */
527 /* ---------------------------------------------------------------------- */
528
529 /* (1) Five calls to UMF_malloc are made, for a total space of
530 * 2 * (n_row + n_col) + 4 integers + sizeof (SymbolicType).
531 * sizeof (SymbolicType) is a small constant. This space is part of the
532 * Symbolic object and is not freed unless an error occurs. If A is square
533 * then this is about 4*n integers.
534 */
535
536 Symbolic = (SymbolicType *) UMF_malloc (1, sizeof (SymbolicType)) ;
537
538 if (!Symbolic)
539 {
540 /* If we fail here, Symbolic is NULL and thus it won't be */
541 /* dereferenced by UMFPACK_free_symbolic, as called by error ( ). */
542 DEBUGm4 (("out of memory: symbolic object\n")) ;
543 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
544 error (&Symbolic, (SWType *) NULL) ;
545 return (UMFPACK_ERROR_out_of_memory) ;
546 }
547
548 /* We now know that Symbolic has been allocated */
549 Symbolic->valid = 0 ;
550 Symbolic->Chain_start = (Int *) NULL ;
551 Symbolic->Chain_maxrows = (Int *) NULL ;
552 Symbolic->Chain_maxcols = (Int *) NULL ;
553 Symbolic->Front_npivcol = (Int *) NULL ;
554 Symbolic->Front_parent = (Int *) NULL ;
555 Symbolic->Front_1strow = (Int *) NULL ;
556 Symbolic->Front_leftmostdesc = (Int *) NULL ;
557 Symbolic->Esize = (Int *) NULL ;
558 Symbolic->esize = 0 ;
559
560 Symbolic->Cperm_init = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
561 Symbolic->Rperm_init = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
562 Symbolic->Cdeg = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
563 Symbolic->Rdeg = (Int *) UMF_malloc (n_row+1, sizeof (Int)) ;
564 Symbolic->Diagonal_map = (Int *) NULL ;
565
566 Cperm_init = Symbolic->Cperm_init ;
567 Rperm_init = Symbolic->Rperm_init ;
568 Cdeg = Symbolic->Cdeg ;
569 Rdeg = Symbolic->Rdeg ;
570
571 if (!Cperm_init || !Rperm_init || !Cdeg || !Rdeg)
572 {
573 DEBUGm4 (("out of memory: symbolic perm\n")) ;
574 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
575 error (&Symbolic, (SWType *) NULL) ;
576 return (UMFPACK_ERROR_out_of_memory) ;
577 }
578
579 Symbolic->n_row = n_row ;
580 Symbolic->n_col = n_col ;
581 Symbolic->nz = nz ;
582 Symbolic->nb = nb ;
583
584 /* ---------------------------------------------------------------------- */
585 /* check user's input permutation */
586 /* ---------------------------------------------------------------------- */
587
588 if (Quser != (Int *) NULL)
589 {
590 /* use Cperm_init as workspace to check input permutation */
591 if (!UMF_is_permutation (Quser, Cperm_init, n_col, n_col))
592 {
593 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_permutation ;
594 error (&Symbolic, (SWType *) NULL) ;
595 return (UMFPACK_ERROR_invalid_permutation) ;
596 }
597 }
598
599 /* ---------------------------------------------------------------------- */
600 /* allocate workspace */
601 /* ---------------------------------------------------------------------- */
602
603 /* (2) Eleven calls to UMF_malloc are made, for workspace of size
604 * Clen + nz + 7*n_col + 2*n_row + 2 integers. Clen is the larger of
605 * MAX (2*nz, 4*n_col) + 8*n_col + 6*n_row + n_col + nz/5 and
606 * 2.4*nz + 8 * MIN (n_row, n_col) + MAX (n_row, n_col, nz)
607 * If A is square and non-singular, then Clen is
608 * MAX (MAX (2*nz, 4*n) + 7*n + nz/5, 3.4*nz) + 8*n
609 * If A has at least 4*n nonzeros then Clen is
610 * MAX (2.2*nz + 7*n, 3.4*nz) + 8*n
611 * If A has at least (7/1.2)*n nonzeros, (about 5.8*n), then Clen is
612 * 3.4*nz + 8*n
613 * This space will be free'd when this routine finishes.
614 *
615 * Total space thus far is about 3.4nz + 12n integers.
616 * For the double precision, 32-bit integer version, the user's matrix
617 * requires an equivalent space of 3*nz + n integers. So this space is just
618 * slightly larger than the user's input matrix (including the numerical
619 * values themselves).
620 */
621
622 SW = &SWspace ; /* used for UMFPACK_symbolic only */
623
624 /* Note that SW->Front_* does not include the dummy placeholder front. */
625 /* This space is accounted for by the SYM_WORK_USAGE macro. */
626
627 /* this is free'd early */
628 SW->Si = (Int *) UMF_malloc (nz, sizeof (Int)) ;
629 SW->Sp = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
630 SW->InvRperm1 = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
631 SW->Cperm1 = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
632
633 /* this is free'd late */
634 SW->Ci = (Int *) UMF_malloc (Clen, sizeof (Int)) ;
635 SW->Front_npivcol = (Int *) UMF_malloc (n_col + 1, sizeof (Int)) ;
636 SW->Front_nrows = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
637 SW->Front_ncols = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
638 SW->Front_parent = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
639 SW->Front_cols = (Int *) UMF_malloc (n_col, sizeof (Int)) ;
640 SW->Rperm1 = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
641 SW->InFront = (Int *) UMF_malloc (n_row, sizeof (Int)) ;
642
643 /* this is allocated later, and free'd after Cperm1 but before Ci */
644 SW->Rperm_2by2 = (Int *) NULL ; /* will be nn Int's */
645
646 /* this is allocated last, and free'd first */
647 SW->Rs = (double *) NULL ; /* will be n_row double's */
648
649 Ci = SW->Ci ;
650 Fr_npivcol = SW->Front_npivcol ;
651 Fr_nrows = SW->Front_nrows ;
652 Fr_ncols = SW->Front_ncols ;
653 Fr_parent = SW->Front_parent ;
654 Fr_cols = SW->Front_cols ;
655 Cperm1 = SW->Cperm1 ;
656 Rperm1 = SW->Rperm1 ;
657 Si = SW->Si ;
658 Sp = SW->Sp ;
659 InvRperm1 = SW->InvRperm1 ;
660 Rperm_2by2 = (Int *) NULL ;
661 InFront = SW->InFront ;
662
663 if (!Ci || !Fr_npivcol || !Fr_nrows || !Fr_ncols || !Fr_parent || !Fr_cols
664 || !Cperm1 || !Rperm1 || !Si || !Sp || !InvRperm1 || !InFront)
665 {
666 DEBUGm4 (("out of memory: symbolic work\n")) ;
667 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
668 error (&Symbolic, SW) ;
669 return (UMFPACK_ERROR_out_of_memory) ;
670 }
671
672 DEBUG0 (("Symbolic UMF_malloc_count - init_count = " ID "\n",
673 UMF_malloc_count - init_count)) ;
674 ASSERT (UMF_malloc_count == init_count + 17) ;
675
676 /* ---------------------------------------------------------------------- */
677 /* find the row and column singletons */
678 /* ---------------------------------------------------------------------- */
679
680 /* [ use first nz + n_row + MAX (n_row, n_col) entries in Ci as workspace,
681 * and use Rperm_init as workspace */
682 ASSERT (Clen >= nz + n_row + MAX (n_row, n_col)) ;
683
684 status = UMF_singletons (n_row, n_col, Ap, Ai, Quser, strategy,
685 Cdeg, Cperm1, Rdeg,
686 Rperm1, InvRperm1, &n1, &n1c, &n1r, &nempty_col, &nempty_row, &is_sym,
687 &max_rdeg, /* workspace: */ Rperm_init, Ci, Ci + nz, Ci + nz + n_row) ;
688
689 /* ] done using Rperm_init and Ci as workspace */
690
691 /* InvRperm1 is now the inverse of Rperm1 */
692
693 if (status != UMFPACK_OK)
694 {
695 DEBUGm4 (("matrix invalid: UMF_singletons\n")) ;
696 Info [UMFPACK_STATUS] = status ;
697 error (&Symbolic, SW) ;
698 return (status) ;
699 }
700 Info [UMFPACK_NEMPTY_COL] = nempty_col ;
701 Info [UMFPACK_NEMPTY_ROW] = nempty_row ;
702 Info [UMFPACK_NDENSE_COL] = 0 ; /* # dense rows/cols recomputed below */
703 Info [UMFPACK_NDENSE_ROW] = 0 ;
704 Info [UMFPACK_COL_SINGLETONS] = n1c ;
705 Info [UMFPACK_ROW_SINGLETONS] = n1r ;
706 Info [UMFPACK_S_SYMMETRIC] = is_sym ;
707
708 nempty = MIN (nempty_col, nempty_row) ;
709 Symbolic->nempty_row = nempty_row ;
710 Symbolic->nempty_col = nempty_col ;
711
712 /* UMF_singletons has verified that the user's input matrix is valid */
713 ASSERT (AMD_valid (n_row, n_col, Ap, Ai) == AMD_OK) ;
714
715 Symbolic->n1 = n1 ;
716 Symbolic->nempty = nempty ;
717 ASSERT (n1 <= n_inner) ;
718 n2 = nn - n1 - nempty ;
719
720 dense_row_threshold =
721 UMFPACK_DENSE_DEGREE_THRESHOLD (drow, n_col - n1 - nempty_col) ;
722 Symbolic->dense_row_threshold = dense_row_threshold ;
723
724 if (!is_sym)
725 {
726 /* either the pruned submatrix rectangular, or it is square and
727 * Rperm [n1 .. n-nempty-1] is not the same as Cperm [n1 .. n-nempty-1].
728 * Switch to the unsymmetric strategy, ignoring user-requested
729 * strategy. */
730 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
731 DEBUGm4 (("Strategy: Unsymmetric singletons\n")) ;
732 }
733
734 /* ---------------------------------------------------------------------- */
735 /* determine symmetry, nzdiag, and degrees of S+S' */
736 /* ---------------------------------------------------------------------- */
737
738 /* S is the matrix obtained after removing singletons
739 * = A (Cperm1 [n1..n_col-nempty_col-1], Rperm1 [n1..n_row-nempty_row-1])
740 */
741
742 Wq = Rperm_init ; /* use Rperm_init as workspace for Wq [ */
743 Sdeg = Cperm_init ; /* use Cperm_init as workspace for Sdeg [ */
744 sym = EMPTY ;
745 nzaat = EMPTY ;
746 nzdiag = EMPTY ;
747 for (i = 0 ; i < AMD_INFO ; i++)
748 {
749 amd_Info [i] = EMPTY ;
750 }
751
752 if (strategy != UMFPACK_STRATEGY_UNSYMMETRIC)
753 {
754 /* This also determines the degree of each node in S+S' (Sdeg), which
755 * is needed by the 2-by-2 strategy, the symmetry of S, and the number
756 * of nonzeros on the diagonal of S. */
757 ASSERT (n_row == n_col) ;
758 ASSERT (nempty_row == nempty_col) ;
759
760 /* get the count of nonzeros on the diagonal of S, excluding explicitly
761 * zero entries. nzdiag = amd_Info [AMD_NZDIAG] counts the zero entries
762 * in S. */
763
764 nzdiag = prune_singletons (n1, nn, Ap, Ai, Ax,
765 #ifdef COMPLEX
766 Az,
767 #endif
768 Cperm1, InvRperm1, Si, Sp
769 #ifndef NDEBUG
770 , Rperm1, nn
771 #endif
772 ) ;
773
774 /* use Ci as workspace to sort S into R, if needed [ */
775 if (Quser != (Int *) NULL)
776 {
777 /* need to sort the columns of S first */
778 Rp = Ci ;
779 Ri = Ci + (n_row) + 1 ;
780 (void) UMF_transpose (n2, n2, Sp, Si, (double *) NULL,
781 (Int *) NULL, (Int *) NULL, 0,
782 Rp, Ri, (double *) NULL, Wq, FALSE
783 #ifdef COMPLEX
784 , (double *) NULL, (double *) NULL, FALSE
785 #endif
786 ) ;
787 }
788 else
789 {
790 /* S already has sorted columns */
791 Rp = Sp ;
792 Ri = Si ;
793 }
794 ASSERT (AMD_valid (n2, n2, Rp, Ri) == AMD_OK) ;
795
796 nzaat = AMD_aat (n2, Rp, Ri, Sdeg, Wq, amd_Info) ;
797 sym = amd_Info [AMD_SYMMETRY] ;
798 Info [UMFPACK_N2] = n2 ;
799 /* nzdiag = amd_Info [AMD_NZDIAG] counts the zero entries of S too */
800
801 /* done using Ci as workspace to sort S into R ] */
802
803 #ifndef NDEBUG
804 for (k = 0 ; k < n2 ; k++)
805 {
806 ASSERT (Sdeg [k] >= 0 && Sdeg [k] < n2) ;
807 }
808 ASSERT (Sp [n2] - n2 <= nzaat && nzaat <= 2 * Sp [n2]) ;
809 DEBUG0 (("Explicit zeros: " ID " %g\n", nzdiag, amd_Info [AMD_NZDIAG])) ;
810 #endif
811
812 }
813
814 /* get statistics from amd_aat, if computed */
815 Symbolic->sym = sym ;
816 Symbolic->nzaat = nzaat ;
817 Symbolic->nzdiag = nzdiag ;
818 Symbolic->amd_dmax = EMPTY ;
819
820 Info [UMFPACK_PATTERN_SYMMETRY] = sym ;
821 Info [UMFPACK_NZ_A_PLUS_AT] = nzaat ;
822 Info [UMFPACK_NZDIAG] = nzdiag ;
823
824 /* ---------------------------------------------------------------------- */
825 /* determine the initial strategy based on symmetry and nnz (diag (S)) */
826 /* ---------------------------------------------------------------------- */
827
828 if (strategy == UMFPACK_STRATEGY_AUTO)
829 {
830 if (sym < 0.10)
831 {
832 /* highly unsymmetric: use the unsymmetric strategy */
833 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
834 DEBUGm4 (("Strategy: select unsymmetric\n")) ;
835 }
836 else if (sym >= 0.7 && nzdiag == n2)
837 {
838 /* mostly symmetric, zero-free diagonal: use symmetric strategy */
839 strategy = UMFPACK_STRATEGY_SYMMETRIC ;
840 DEBUGm4 (("Strategy: select symmetric\n")) ;
841 }
842 else
843 {
844 /* Evaluate the symmetric 2-by-2 strategy, and select it, or
845 * the unsymmetric strategy if the 2-by-2 strategy doesn't look
846 * promising. */
847 strategy = UMFPACK_STRATEGY_2BY2 ;
848 DEBUGm4 (("Strategy: try 2-by-2\n")) ;
849 }
850 }
851
852 /* ---------------------------------------------------------------------- */
853 /* try the 2-by-2 strategy */
854 /* ---------------------------------------------------------------------- */
855
856 /* (3) If the 2-by-2 strategy is attempted, additional workspace of size
857 * nn integers and nn double's is allocated, where nn = n_row = n_col.
858 * The real workspace is immediately free'd. The integer workspace of
859 * size nn remains until the end of umfpack_qsymbolic. */
860
861 /* If the resulting matrix S (Rperm_2by2, :) is too unsymmetric, then the
862 * unsymmetric strategy will be used instead. */
863
864 if (strategy == UMFPACK_STRATEGY_2BY2)
865 {
866 double sym2 ;
867 Int *Blen, *W, nz_papat, nzd2, nweak, unmatched, Clen3 ;
868
869 /* ------------------------------------------------------------------ */
870 /* get workspace for UMF_2by2 */
871 /* ------------------------------------------------------------------ */
872
873 ASSERT (n_row == n_col && nn == n_row) ;
874
875 #ifndef NDEBUG
876 for (k = 0 ; k < n2 ; k++)
877 {
878 ASSERT (Sdeg [k] >= 0 && Sdeg [k] < n2) ;
879 }
880 #endif
881
882 /* allocate Rperm_2by2 */
883 SW->Rperm_2by2 = (Int *) UMF_malloc (nn, sizeof (Int)) ;
884 Rperm_2by2 = SW->Rperm_2by2 ;
885 if (Rperm_2by2 == (Int *) NULL)
886 {
887 DEBUGm4 (("out of memory: Rperm_2by2\n")) ;
888 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
889 error (&Symbolic, SW) ;
890 return (UMFPACK_ERROR_out_of_memory) ;
891 }
892
893 /* allocate Ri from the tail end of Ci [ */
894 Clen3 = Clen - (MAX (nn, nz) + 1) ;
895 Ri = Ci + Clen3 ;
896 ASSERT (Clen3 >= nz) ; /* space required for UMF_2by2 */
897
898 /* use Fr_* as workspace for Rp, Blen, and W [ */
899 Rp = Fr_npivcol ;
900 Blen = Fr_ncols ;
901 W = Fr_cols ;
902
903 if (scale != UMFPACK_SCALE_NONE)
904 {
905 SW->Rs = (double *) UMF_malloc (nn, sizeof (double)) ;
906 if (SW->Rs == (double *) NULL)
907 {
908 DEBUGm4 (("out of memory: scale factors for 2-by-2\n")) ;
909 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
910 error (&Symbolic, SW) ;
911 return (UMFPACK_ERROR_out_of_memory) ;
912 }
913 }
914
915 /* ------------------------------------------------------------------ */
916 /* find the 2-by-2 row permutation */
917 /* ------------------------------------------------------------------ */
918
919 /* find a row permutation Rperm_2by2 such that S (Rperm_2by2, :)
920 * has a healthy diagonal */
921
922 UMF_2by2 (nn, Ap, Ai, Ax,
923 #ifdef COMPLEX
924 Az,
925 #endif
926 tol, scale, Cperm1,
927 #ifndef NDEBUG
928 Rperm1,
929 #endif
930 InvRperm1, n1, nempty, Sdeg, Rperm_2by2, &nweak, &unmatched,
931 Ri, Rp, SW->Rs, Blen, W, Ci, Wq) ;
932 DEBUGm3 (("2by2: nweak " ID " unmatched " ID "\n", nweak, unmatched)) ;
933 Info [UMFPACK_2BY2_NWEAK] = nweak ;
934 Info [UMFPACK_2BY2_UNMATCHED] = unmatched ;
935
936 SW->Rs = (double *) UMF_free ((void *) SW->Rs) ;
937
938 /* R = S (Rperm_2by2,:)' */
939 (void) UMF_transpose (n2, n2, Sp, Si, (double *) NULL, Rperm_2by2,
940 (Int *) NULL, 0, Rp, Ri, (double *) NULL, W, FALSE
941 #ifdef COMPLEX
942 , (double *) NULL, (double *) NULL, FALSE
943 #endif
944 ) ;
945 ASSERT (AMD_valid (n2, n2, Rp, Ri) == AMD_OK) ;
946
947 /* contents of Si and Sp no longer needed, but the space is
948 * still needed */
949
950 /* ------------------------------------------------------------------ */
951 /* find symmetry of S (Rperm_2by2, :)', and prepare to order with AMD */
952 /* ------------------------------------------------------------------ */
953
954 for (i = 0 ; i < AMD_INFO ; i++)
955 {
956 amd_Info [i] = EMPTY ;
957 }
958 nz_papat = AMD_aat (n2, Rp, Ri, Sdeg, Wq, amd_Info) ;
959 sym2 = amd_Info [AMD_SYMMETRY] ;
960 nzd2 = amd_Info [AMD_NZDIAG] ;
961
962 Info [UMFPACK_2BY2_PATTERN_SYMMETRY] = sym2 ;
963 Info [UMFPACK_2BY2_NZ_PA_PLUS_PAT] = nz_papat ;
964 Info [UMFPACK_2BY2_NZDIAG] = nzd2 ;
965
966 DEBUG0 (("2by2: sym2 %g nzd2 " ID " n2 " ID "\n", sym2, nzd2, n2)) ;
967
968 /* ------------------------------------------------------------------ */
969 /* evaluate the 2-by-2 results */
970 /* ------------------------------------------------------------------ */
971
972 if (user_auto_strategy)
973 {
974 if ((sym2 > 1.1 * sym) && (nzd2 > 0.9 * n2))
975 {
976 /* 2-by-2 made it much more symmetric */
977 DEBUGm4 (("eval Strategy 2by2: much more symmetric: 2by2\n")) ;
978 strategy = UMFPACK_STRATEGY_2BY2 ;
979 }
980 else if (sym2 < 0.7 * sym)
981 {
982 /* 2-by-2 made it much more unsymmetric */
983 DEBUGm4 (("eval Strategy 2by2: much more UNsymmetric:unsym\n"));
984 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
985 }
986 else if (sym2 < 0.25)
987 {
988 DEBUGm4 (("eval Strategy 2by2: is UNsymmetric: unsym\n"));
989 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
990 }
991 else if (sym2 >= 0.51)
992 {
993 DEBUGm4 (("eval Strategy 2by2: sym2 >= 0.51: 2by2\n")) ;
994 strategy = UMFPACK_STRATEGY_2BY2 ;
995 }
996 else if (sym2 >= 0.999 * sym)
997 {
998 /* 2-by-2 improved symmetry, or made it only slightly worse */
999 DEBUGm4 (("eval Strategy 2by2: sym2 >= 0.999 sym: 2by2\n")) ;
1000 strategy = UMFPACK_STRATEGY_2BY2 ;
1001 }
1002 else
1003 {
1004 /* can't decide what to do, so pick the unsymmetric strategy */
1005 DEBUGm4 (("eval Strategy 2by2: punt: unsym\n"));
1006 strategy = UMFPACK_STRATEGY_UNSYMMETRIC ;
1007 }
1008 }
1009
1010 /* ------------------------------------------------------------------ */
1011 /* if the 2-by-2 strategy is selected: */
1012 /* ------------------------------------------------------------------ */
1013
1014 if (strategy == UMFPACK_STRATEGY_2BY2)
1015 {
1016 if (Quser == (Int *) NULL)
1017 {
1018 /* 2-by-2 strategy is successful */
1019 /* compute amd (S) */
1020 Int *Qinv = Fr_npivcol ;
1021 ASSERT (Clen3 >= (nz_papat + nz_papat/5 + nn) + 7*nn) ;
1022 do_amd (n2, Rp, Ri, Wq, Qinv, Sdeg, Clen3, Ci,
1023 amd_Control, amd_Info, Symbolic, Info) ;
1024 /* combine the singleton ordering and the AMD ordering */
1025 combine_ordering (n1, nempty, nn, Cperm_init, Cperm1, Qinv) ;
1026 }
1027 /* fix Rperm_2by2 to reflect A, not S */
1028 for (k = 0 ; k < n1 ; k++)
1029 {
1030 oldcol = Cperm1 [k] ;
1031 i = k ;
1032 oldrow = Rperm1 [k] ;
1033 W [oldcol] = oldrow ;
1034 }
1035 for (k = n1 ; k < nn - nempty ; k++)
1036 {
1037 oldcol = Cperm1 [k] ;
1038 i = Rperm_2by2 [k - n1] + n1 ;
1039 oldrow = Rperm1 [i] ;
1040 W [oldcol] = oldrow ;
1041 }
1042 for (k = nn - nempty ; k < nn ; k++)
1043 {
1044 oldcol = Cperm1 [k] ;
1045 i = k ;
1046 oldrow = Rperm1 [k] ;
1047 W [oldcol] = oldrow ;
1048 }
1049 for (k = 0 ; k < nn ; k++)
1050 {
1051 Rperm_2by2 [k] = W [k] ;
1052 }
1053
1054 /* Now, the "diagonal" entry in oldcol (where oldcol is the user's
1055 * name for a column, is the entry in row oldrow (where oldrow is
1056 * the user's name for a row, and oldrow = Rperm_2by2 [oldcol] */
1057 }
1058
1059 /* Fr_* no longer needed for Rp, Blen, W ] */
1060 }
1061
1062 /* ---------------------------------------------------------------------- */
1063 /* finalize the strategy, including fixQ and prefer_diagonal */
1064 /* ---------------------------------------------------------------------- */
1065
1066 if (strategy == UMFPACK_STRATEGY_SYMMETRIC)
1067 {
1068 /* use given Quser or AMD (A+A'), fix Q during factorization,
1069 * prefer diagonal */
1070 DEBUG0 (("\nStrategy: symmetric\n")) ;
1071 ASSERT (n_row == n_col) ;
1072 Symbolic->ordering = UMFPACK_ORDERING_AMD ;
1073 fixQ = TRUE ;
1074 prefer_diagonal = TRUE ;
1075 }
1076 else if (strategy == UMFPACK_STRATEGY_2BY2)
1077 {
1078 /* use Q = given Quser or Q = AMD (PA+PA'), fix Q during factorization,
1079 * prefer diagonal, and factorize PAQ, where P is found by UMF_2by2. */
1080 DEBUG0 (("\nStrategy: symmetric 2-by-2\n")) ;
1081 ASSERT (n_row == n_col) ;
1082 Symbolic->ordering = UMFPACK_ORDERING_AMD ;
1083 fixQ = TRUE ;
1084 prefer_diagonal = TRUE ;
1085 }
1086 else
1087 {
1088 /* use given Quser or COLAMD (A), refine Q during factorization,
1089 * no diagonal preference */
1090 ASSERT (strategy == UMFPACK_STRATEGY_UNSYMMETRIC) ;
1091 DEBUG0 (("\nStrategy: unsymmetric\n")) ;
1092 Symbolic->ordering = UMFPACK_ORDERING_COLAMD ;
1093 fixQ = FALSE ;
1094 prefer_diagonal = FALSE ;
1095 }
1096
1097 if (Quser != (Int *) NULL)
1098 {
1099 Symbolic->ordering = UMFPACK_ORDERING_GIVEN ;
1100 }
1101
1102 if (force_fixQ > 0)
1103 {
1104 fixQ = TRUE ;
1105 DEBUG0 (("Force fixQ true\n")) ;
1106 }
1107 else if (force_fixQ < 0)
1108 {
1109 fixQ = FALSE ;
1110 DEBUG0 (("Force fixQ false\n")) ;
1111 }
1112
1113 DEBUG0 (("Strategy: ordering: " ID "\n", Symbolic->ordering)) ;
1114 DEBUG0 (("Strategy: fixQ: " ID "\n", fixQ)) ;
1115 DEBUG0 (("Strategy: prefer diag " ID "\n", prefer_diagonal)) ;
1116
1117 /* get statistics from amd_aat, if computed */
1118 Symbolic->strategy = strategy ;
1119 Symbolic->fixQ = fixQ ;
1120 Symbolic->prefer_diagonal = prefer_diagonal ;
1121
1122 Info [UMFPACK_STRATEGY_USED] = strategy ;
1123 Info [UMFPACK_ORDERING_USED] = Symbolic->ordering ;
1124 Info [UMFPACK_QFIXED] = fixQ ;
1125 Info [UMFPACK_DIAG_PREFERRED] = prefer_diagonal ;
1126
1127 /* ---------------------------------------------------------------------- */
1128 /* get the AMD ordering for the symmetric strategy */
1129 /* ---------------------------------------------------------------------- */
1130
1131 if (strategy == UMFPACK_STRATEGY_SYMMETRIC && Quser == (Int *) NULL)
1132 {
1133 /* symmetric strategy for a matrix with mostly symmetric pattern */
1134 Int *Qinv = Fr_npivcol ;
1135 ASSERT (n_row == n_col && nn == n_row) ;
1136 ASSERT (Clen >= (nzaat + nzaat/5 + nn) + 7*nn) ;
1137 do_amd (n2, Sp, Si, Wq, Qinv, Sdeg, Clen, Ci,
1138 amd_Control, amd_Info, Symbolic, Info) ;
1139 /* combine the singleton ordering and the AMD ordering */
1140 combine_ordering (n1, nempty, nn, Cperm_init, Cperm1, Qinv) ;
1141 }
1142 /* Sdeg no longer needed ] */
1143 /* done using Rperm_init as workspace for Wq ] */
1144
1145 /* Contents of Si and Sp no longer needed, but the space is still needed */
1146
1147 /* ---------------------------------------------------------------------- */
1148 /* use the user's input column ordering (already in Cperm1) */
1149 /* ---------------------------------------------------------------------- */
1150
1151 if (Quser != (Int *) NULL)
1152 {
1153 for (k = 0 ; k < n_col ; k++)
1154 {
1155 Cperm_init [k] = Cperm1 [k] ;
1156 }
1157 }
1158
1159 /* ---------------------------------------------------------------------- */
1160 /* use COLAMD to order the matrix */
1161 /* ---------------------------------------------------------------------- */
1162
1163 if (strategy == UMFPACK_STRATEGY_UNSYMMETRIC && Quser == (Int *) NULL)
1164 {
1165
1166 /* ------------------------------------------------------------------ */
1167 /* copy the matrix into colamd workspace (colamd destroys its input) */
1168 /* ------------------------------------------------------------------ */
1169
1170 /* C = A (Cperm1 (n1+1:end), Rperm1 (n1+1:end)), where Ci is used as
1171 * the row indices and Cperm_init (on input) is used as the column
1172 * pointers. */
1173
1174 (void) prune_singletons (n1, n_col, Ap, Ai,
1175 (double *) NULL,
1176 #ifdef COMPLEX
1177 (double *) NULL,
1178 #endif
1179 Cperm1, InvRperm1, Ci, Cperm_init
1180 #ifndef NDEBUG
1181 , Rperm1, n_row
1182 #endif
1183 ) ;
1184
1185 /* ------------------------------------------------------------------ */
1186 /* set UMF_colamd defaults */
1187 /* ------------------------------------------------------------------ */
1188
1189 UMF_colamd_set_defaults (knobs) ;
1190 knobs [COLAMD_DENSE_ROW] = drow ;
1191 knobs [COLAMD_DENSE_COL] = dcol ;
1192 knobs [COLAMD_AGGRESSIVE] = aggressive ;
1193
1194 /* ------------------------------------------------------------------ */
1195 /* check input matrix and find the initial column pre-ordering */
1196 /* ------------------------------------------------------------------ */
1197
1198 /* NOTE: umf_colamd is not given any original empty rows or columns.
1199 * Those have already been removed via prune_singletons, above. The
1200 * umf_colamd routine has been modified to assume that all rows and
1201 * columns have at least one entry in them. It will break if it is
1202 * given empty rows or columns (an assertion is triggered when running
1203 * in debug mode. */
1204
1205 (void) UMF_colamd (
1206 n_row - n1 - nempty_row,
1207 n_col - n1 - nempty_col,
1208 Clen, Ci, Cperm_init, knobs, colamd_stats,
1209 Fr_npivcol, Fr_nrows, Fr_ncols, Fr_parent, Fr_cols, &nfr,
1210 InFront) ;
1211 ASSERT (colamd_stats [COLAMD_EMPTY_ROW] == 0) ;
1212 ASSERT (colamd_stats [COLAMD_EMPTY_COL] == 0) ;
1213
1214 /* # of dense rows will be recomputed below */
1215 Info [UMFPACK_NDENSE_ROW] = colamd_stats [COLAMD_DENSE_ROW] ;
1216 Info [UMFPACK_NDENSE_COL] = colamd_stats [COLAMD_DENSE_COL] ;
1217 Info [UMFPACK_SYMBOLIC_DEFRAG] = colamd_stats [COLAMD_DEFRAG_COUNT] ;
1218
1219 /* re-analyze if any "dense" rows or cols ignored by UMF_colamd */
1220 do_UMF_analyze =
1221 colamd_stats [COLAMD_DENSE_ROW] > 0 ||
1222 colamd_stats [COLAMD_DENSE_COL] > 0 ;
1223
1224 /* Combine the singleton and colamd ordering into Cperm_init */
1225 /* Note that colamd returns its inverse permutation in Ci */
1226 combine_ordering (n1, nempty_col, n_col, Cperm_init, Cperm1, Ci) ;
1227
1228 /* contents of Ci no longer needed */
1229
1230 #ifndef NDEBUG
1231 for (col = 0 ; col < n_col ; col++)
1232 {
1233 DEBUG1 (("Cperm_init [" ID "] = " ID "\n", col, Cperm_init[col]));
1234 }
1235 /* make sure colamd returned a valid permutation */
1236 ASSERT (Cperm_init != (Int *) NULL) ;
1237 ASSERT (UMF_is_permutation (Cperm_init, Ci, n_col, n_col)) ;
1238 #endif
1239
1240 }
1241 else
1242 {
1243
1244 /* ------------------------------------------------------------------ */
1245 /* do not call colamd - use input Quser or AMD instead */
1246 /* ------------------------------------------------------------------ */
1247
1248 /* The ordering (Quser or Qamd) is already in Cperm_init */
1249 do_UMF_analyze = TRUE ;
1250
1251 }
1252
1253 Cperm_init [n_col] = EMPTY ; /* unused in Cperm_init */
1254
1255 /* ---------------------------------------------------------------------- */
1256 /* AMD ordering, if it exists, has been copied into Cperm_init */
1257 /* ---------------------------------------------------------------------- */
1258
1259 #ifndef NDEBUG
1260 DEBUG3 (("Cperm_init column permutation:\n")) ;
1261 ASSERT (UMF_is_permutation (Cperm_init, Ci, n_col, n_col)) ;
1262 for (k = 0 ; k < n_col ; k++)
1263 {
1264 DEBUG3 ((ID "\n", Cperm_init [k])) ;
1265 }
1266 /* ensure that empty columns have been placed last in A (:,Cperm_init) */
1267 for (newj = 0 ; newj < n_col ; newj++)
1268 {
1269 /* empty columns will be last in A (:, Cperm_init (1:n_col)) */
1270 j = Cperm_init [newj] ;
1271 ASSERT (IMPLIES (newj >= n_col-nempty_col, Cdeg [j] == 0)) ;
1272 ASSERT (IMPLIES (newj < n_col-nempty_col, Cdeg [j] > 0)) ;
1273 }
1274 #endif
1275
1276 /* ---------------------------------------------------------------------- */
1277 /* symbolic factorization (unless colamd has already done it) */
1278 /* ---------------------------------------------------------------------- */
1279
1280 if (do_UMF_analyze)
1281 {
1282
1283 Int *W, *Bp, *Bi, *Cperm2, ok, *P, Clen2, bsize, Clen0 ;
1284
1285 /* ------------------------------------------------------------------ */
1286 /* construct column pre-ordered, pruned submatrix */
1287 /* ------------------------------------------------------------------ */
1288
1289 /* S = column form submatrix after removing singletons and applying
1290 * initial column ordering (includes singleton ordering) */
1291 (void) prune_singletons (n1, n_col, Ap, Ai,
1292 (double *) NULL,
1293 #ifdef COMPLEX
1294 (double *) NULL,
1295 #endif
1296 Cperm_init, InvRperm1, Si, Sp
1297 #ifndef NDEBUG
1298 , Rperm1, n_row
1299 #endif
1300 ) ;
1301
1302 /* ------------------------------------------------------------------ */
1303 /* Ci [0 .. Clen-1] holds the following work arrays:
1304
1305 first Clen0 entries empty space, where Clen0 =
1306 Clen - (nn+1 + 2*nn + n_col)
1307 and Clen0 >= nz + n_col
1308 next nn+1 entries Bp [0..nn]
1309 next nn entries Link [0..nn-1]
1310 next nn entries W [0..nn-1]
1311 last n_col entries Cperm2 [0..n_col-1]
1312
1313 We have Clen >= n_col + MAX (nz,n_col) + 3*nn+1 + n_col,
1314 So Clen0 >= 2*n_col as required for AMD_postorder
1315 and Clen0 >= n_col + nz as required
1316 */
1317
1318 Clen0 = Clen - (nn+1 + 2*nn + n_col) ;
1319 Bp = Ci + Clen0 ;
1320 Link = Bp + (nn+1) ;
1321 W = Link + nn ;
1322 Cperm2 = W + nn ;
1323 ASSERT (Cperm2 + n_col == Ci + Clen) ;
1324 ASSERT (Clen0 >= nz + n_col) ;
1325 ASSERT (Clen0 >= 2*n_col) ;
1326
1327 /* ------------------------------------------------------------------ */
1328 /* P = order that rows will be used in UMF_analyze */
1329 /* ------------------------------------------------------------------ */
1330
1331 /* use W to mark rows, and use Link for row permutation P [ [ */
1332 for (row = 0 ; row < n_row - n1 ; row++)
1333 {
1334 W [row] = FALSE ;
1335 }
1336 P = Link ;
1337
1338 k = 0 ;
1339
1340 for (col = 0 ; col < n_col - n1 ; col++)
1341 {
1342 /* empty columns are last in S */
1343 for (p = Sp [col] ; p < Sp [col+1] ; p++)
1344 {
1345 row = Si [p] ;
1346 if (!W [row])
1347 {
1348 /* this row has just been seen for the first time */
1349 W [row] = TRUE ;
1350 P [k++] = row ;
1351 }
1352 }
1353 }
1354
1355 /* If the matrix has truly empty rows, then P will not be */
1356 /* complete, and visa versa. The matrix is structurally singular. */
1357 nempty_row = n_row - n1 - k ;
1358 if (k < n_row - n1)
1359 {
1360 /* complete P by putting empty rows last in their natural order, */
1361 /* rather than declaring an error (the matrix is singular) */
1362 for (row = 0 ; row < n_row - n1 ; row++)
1363 {
1364 if (!W [row])
1365 {
1366 /* W [row] = TRUE ; (not required) */
1367 P [k++] = row ;
1368 }
1369 }
1370 }
1371
1372 /* contents of W no longer needed ] */
1373
1374 #ifndef NDEBUG
1375 DEBUG3 (("Induced row permutation:\n")) ;
1376 ASSERT (k == n_row - n1) ;
1377 ASSERT (UMF_is_permutation (P, W, n_row - n1, n_row - n1)) ;
1378 for (k = 0 ; k < n_row - n1 ; k++)
1379 {
1380 DEBUG3 ((ID "\n", P [k])) ;
1381 }
1382 #endif
1383
1384 /* ------------------------------------------------------------------ */
1385 /* B = row-form of the pattern of S (excluding empty columns) */
1386 /* ------------------------------------------------------------------ */
1387
1388 /* Ci [0 .. Clen-1] holds the following work arrays:
1389
1390 first Clen2 entries empty space, must be at least >= n_col
1391 next max (nz,1) Bi [0..max (nz,1)-1]
1392 next nn+1 entries Bp [0..nn]
1393 next nn entries Link [0..nn-1]
1394 next nn entries W [0..nn-1]
1395 last n_col entries Cperm2 [0..n_col-1]
1396
1397 This memory usage is accounted for by the UMF_ANALYZE_CLEN
1398 macro.
1399 */
1400
1401 Clen2 = Clen0 ;
1402 snz = Sp [n_col - n1] ;
1403 bsize = MAX (snz, 1) ;
1404 Clen2 -= bsize ;
1405 Bi = Ci + Clen2 ;
1406 ASSERT (Clen2 >= n_col) ;
1407
1408 (void) UMF_transpose (n_row - n1, n_col - n1 - nempty_col,
1409 Sp, Si, (double *) NULL,
1410 P, (Int *) NULL, 0, Bp, Bi, (double *) NULL, W, FALSE
1411 #ifdef COMPLEX
1412 , (double *) NULL, (double *) NULL, FALSE
1413 #endif
1414 ) ;
1415
1416 /* contents of Si and Sp no longer needed */
1417
1418 /* contents of P (same as Link) and W not needed */
1419 /* still need Link and W as work arrays, though ] */
1420
1421 ASSERT (Bp [0] == 0) ;
1422 ASSERT (Bp [n_row - n1] == snz) ;
1423
1424 /* increment Bp to point into Ci, not Bi */
1425 for (i = 0 ; i <= n_row - n1 ; i++)
1426 {
1427 Bp [i] += Clen2 ;
1428 }
1429 ASSERT (Bp [0] == Clen0 - bsize) ;
1430 ASSERT (Bp [n_row - n1] <= Clen0) ;
1431
1432 /* Ci [0 .. Clen-1] holds the following work arrays:
1433
1434 first Clen0 entries Ci [0 .. Clen0-1], where the col indices
1435 of B are at the tail end of this part,
1436 and Bp [0] = Clen2 >= n_col. Note
1437 that Clen0 = Clen2 + max (snz,1).
1438 next nn+1 entries Bp [0..nn]
1439 next nn entries Link [0..nn-1]
1440 next nn entries W [0..nn-1]
1441 last n_col entries Cperm2 [0..n_col-1]
1442 */
1443
1444 /* ------------------------------------------------------------------ */
1445 /* analyze */
1446 /* ------------------------------------------------------------------ */
1447
1448 /* only analyze the non-empty, non-singleton part of the matrix */
1449 ok = UMF_analyze (
1450 n_row - n1 - nempty_row,
1451 n_col - n1 - nempty_col,
1452 Ci, Bp, Cperm2, fixQ, W, Link,
1453 Fr_ncols, Fr_nrows, Fr_npivcol,
1454 Fr_parent, &nfr, &analyze_compactions) ;
1455 if (!ok)
1456 {
1457 /* :: internal error in umf_analyze :: */
1458 Info [UMFPACK_STATUS] = UMFPACK_ERROR_internal_error ;
1459 error (&Symbolic, SW) ;
1460 return (UMFPACK_ERROR_internal_error) ;
1461 }
1462 Info [UMFPACK_SYMBOLIC_DEFRAG] += analyze_compactions ;
1463
1464 /* ------------------------------------------------------------------ */
1465 /* combine the input permutation and UMF_analyze's permutation */
1466 /* ------------------------------------------------------------------ */
1467
1468 if (!fixQ)
1469 {
1470 /* Cperm2 is the column etree post-ordering */
1471 ASSERT (UMF_is_permutation (Cperm2, W,
1472 n_col-n1-nempty_col, n_col-n1-nempty_col)) ;
1473
1474 /* Note that the empty columns remain at the end of Cperm_init */
1475 for (k = 0 ; k < n_col - n1 - nempty_col ; k++)
1476 {
1477 W [k] = Cperm_init [n1 + Cperm2 [k]] ;
1478 }
1479
1480 for (k = 0 ; k < n_col - n1 - nempty_col ; k++)
1481 {
1482 Cperm_init [n1 + k] = W [k] ;
1483 }
1484 }
1485
1486 ASSERT (UMF_is_permutation (Cperm_init, W, n_col, n_col)) ;
1487
1488 }
1489
1490 /* ---------------------------------------------------------------------- */
1491 /* free some of the workspace */
1492 /* ---------------------------------------------------------------------- */
1493
1494 /* (4) The real workspace, Rs, of size n_row doubles has already been
1495 * free'd. An additional workspace of size nz + n_col+1 + n_col integers
1496 * is now free'd as well. */
1497
1498 SW->Si = (Int *) UMF_free ((void *) SW->Si) ;
1499 SW->Sp = (Int *) UMF_free ((void *) SW->Sp) ;
1500 SW->Cperm1 = (Int *) UMF_free ((void *) SW->Cperm1) ;
1501 ASSERT (SW->Rs == (double *) NULL) ;
1502
1503 /* ---------------------------------------------------------------------- */
1504 /* determine the size of the Symbolic object */
1505 /* ---------------------------------------------------------------------- */
1506
1507 /* ---------------------------------------------------------------------- */
1508 /* determine the size of the Symbolic object */
1509 /* ---------------------------------------------------------------------- */
1510
1511 nchains = 0 ;
1512 for (i = 0 ; i < nfr ; i++)
1513 {
1514 if (Fr_parent [i] != i+1)
1515 {
1516 nchains++ ;
1517 }
1518 }
1519
1520 Symbolic->nchains = nchains ;
1521 Symbolic->nfr = nfr ;
1522 Symbolic->esize
1523 = (max_rdeg > dense_row_threshold) ? (n_col - n1 - nempty_col) : 0 ;
1524
1525 /* true size of Symbolic object */
1526 Info [UMFPACK_SYMBOLIC_SIZE] = UMF_symbolic_usage (n_row, n_col, nchains,
1527 nfr, Symbolic->esize, prefer_diagonal) ;
1528
1529 /* actual peak memory usage for UMFPACK_symbolic (actual nfr, nchains) */
1530 Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] =
1531 SYM_WORK_USAGE (n_col, n_row, Clen) + Info [UMFPACK_SYMBOLIC_SIZE] ;
1532 Symbolic->peak_sym_usage = Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] ;
1533
1534 DEBUG0 (("Number of fronts: " ID "\n", nfr)) ;
1535
1536 /* ---------------------------------------------------------------------- */
1537 /* allocate the second part of the Symbolic object (Front_*, Chain_*) */
1538 /* ---------------------------------------------------------------------- */
1539
1540 /* (5) UMF_malloc is called 7 or 8 times, for a total space of
1541 * (4*(nfr+1) + 3*(nchains+1) + esize) integers, where nfr is the total
1542 * number of frontal matrices and nchains is the total number of frontal
1543 * matrix chains, and where nchains <= nfr <= n_col. esize is zero if there
1544 * are no dense rows, or n_col-n1-nempty_col otherwise (n1 is the number of
1545 * singletons and nempty_col is the number of empty columns). This space is
1546 * part of the Symbolic object and is not free'd unless an error occurs.
1547 * This is between 7 and about 8n integers when A is square.
1548 */
1549
1550 /* Note that Symbolic->Front_* does include the dummy placeholder front */
1551 Symbolic->Front_npivcol = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1552 Symbolic->Front_parent = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1553 Symbolic->Front_1strow = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1554 Symbolic->Front_leftmostdesc = (Int *) UMF_malloc (nfr+1, sizeof (Int)) ;
1555 Symbolic->Chain_start = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1556 Symbolic->Chain_maxrows = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1557 Symbolic->Chain_maxcols = (Int *) UMF_malloc (nchains+1, sizeof (Int)) ;
1558
1559 fail = (!Symbolic->Front_npivcol || !Symbolic->Front_parent ||
1560 !Symbolic->Front_1strow || !Symbolic->Front_leftmostdesc ||
1561 !Symbolic->Chain_start || !Symbolic->Chain_maxrows ||
1562 !Symbolic->Chain_maxcols) ;
1563
1564 if (Symbolic->esize > 0)
1565 {
1566 Symbolic->Esize = (Int *) UMF_malloc (Symbolic->esize, sizeof (Int)) ;
1567 fail = fail || !Symbolic->Esize ;
1568 }
1569
1570 if (fail)
1571 {
1572 DEBUGm4 (("out of memory: rest of symbolic object\n")) ;
1573 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
1574 error (&Symbolic, SW) ;
1575 return (UMFPACK_ERROR_out_of_memory) ;
1576 }
1577 DEBUG0 (("Symbolic UMF_malloc_count - init_count = " ID "\n",
1578 UMF_malloc_count - init_count)) ;
1579 ASSERT (UMF_malloc_count == init_count + 21
1580 + (SW->Rperm_2by2 != (Int *) NULL)
1581 + (Symbolic->Esize != (Int *) NULL)) ;
1582
1583 Front_npivcol = Symbolic->Front_npivcol ;
1584 Front_parent = Symbolic->Front_parent ;
1585 Front_1strow = Symbolic->Front_1strow ;
1586 Front_leftmostdesc = Symbolic->Front_leftmostdesc ;
1587
1588 Chain_start = Symbolic->Chain_start ;
1589 Chain_maxrows = Symbolic->Chain_maxrows ;
1590 Chain_maxcols = Symbolic->Chain_maxcols ;
1591
1592 Esize = Symbolic->Esize ;
1593
1594 /* ---------------------------------------------------------------------- */
1595 /* assign rows to fronts */
1596 /* ---------------------------------------------------------------------- */
1597
1598 /* find InFront, unless colamd has already computed it */
1599 if (do_UMF_analyze)
1600 {
1601
1602 DEBUGm4 ((">>>>>>>>>Computing Front_1strow from scratch\n")) ;
1603 /* empty rows go to dummy front nfr */
1604 for (row = 0 ; row < n_row ; row++)
1605 {
1606 InFront [row] = nfr ;
1607 }
1608 /* assign the singleton pivot rows to the "empty" front */
1609 for (k = 0 ; k < n1 ; k++)
1610 {
1611 row = Rperm1 [k] ;
1612 InFront [row] = EMPTY ;
1613 }
1614 DEBUG1 (("Front (EMPTY), singleton nrows " ID " ncols " ID "\n", k, k)) ;
1615 newj = n1 ;
1616 for (i = 0 ; i < nfr ; i++)
1617 {
1618 fpivcol = Fr_npivcol [i] ;
1619 f1rows = 0 ;
1620 /* for all pivot columns in front i */
1621 for (kk = 0 ; kk < fpivcol ; kk++, newj++)
1622 {
1623 j = Cperm_init [newj] ;
1624 ASSERT (IMPLIES (newj >= n_col-nempty_col,
1625 Ap [j+1] - Ap [j] == 0));
1626 for (p = Ap [j] ; p < Ap [j+1] ; p++)
1627 {
1628 row = Ai [p] ;
1629 if (InFront [row] == nfr)
1630 {
1631 /* this row belongs to front i */
1632 DEBUG1 ((" Row " ID " in Front " ID "\n", row, i)) ;
1633 InFront [row] = i ;
1634 f1rows++ ;
1635 }
1636 }
1637 }
1638 Front_1strow [i] = f1rows ;
1639 DEBUG1 ((" Front " ID " has 1strows: " ID " pivcols " ID "\n",
1640 i, f1rows, fpivcol)) ;
1641 }
1642
1643 }
1644 else
1645 {
1646
1647 /* COLAMD has already computed InFront, but it is not yet
1648 * InFront [row] = front i, where row is an original row. It is
1649 * InFront [k-n1] = i for k in the range n1 to n_row-nempty_row,
1650 * and where row = Rperm1 [k]. Need to permute InFront. Also compute
1651 * # of original rows assembled into each front.
1652 * [ use Ci as workspace */
1653 DEBUGm4 ((">>>>>>>>>Computing Front_1strow from colamd's InFront\n")) ;
1654 for (i = 0 ; i <= nfr ; i++)
1655 {
1656 Front_1strow [i] = 0 ;
1657 }
1658 /* assign the singleton pivot rows to "empty" front */
1659 for (k = 0 ; k < n1 ; k++)
1660 {
1661 row = Rperm1 [k] ;
1662 Ci [row] = EMPTY ;
1663 }
1664 /* assign the non-empty rows to the front that assembled them */
1665 for ( ; k < n_row - nempty_row ; k++)
1666 {
1667 row = Rperm1 [k] ;
1668 i = InFront [k - n1] ;
1669 ASSERT (i >= EMPTY && i < nfr) ;
1670 if (i != EMPTY)
1671 {
1672 Front_1strow [i]++ ;
1673 }
1674 /* use Ci as permuted version of InFront */
1675 Ci [row] = i ;
1676 }
1677 /* empty rows go to the "dummy" front */
1678 for ( ; k < n_row ; k++)
1679 {
1680 row = Rperm1 [k] ;
1681 Ci [row] = nfr ;
1682 }
1683 /* permute InFront so that InFront [row] = i if the original row is
1684 * in front i */
1685 for (row = 0 ; row < n_row ; row++)
1686 {
1687 InFront [row] = Ci [row] ;
1688 }
1689 /* ] no longer need Ci as workspace */
1690 }
1691
1692 #ifndef NDEBUG
1693 for (row = 0 ; row < n_row ; row++)
1694 {
1695 if (InFront [row] == nfr)
1696 {
1697 DEBUG1 ((" Row " ID " in Dummy Front " ID "\n", row, nfr)) ;
1698 }
1699 else if (InFront [row] == EMPTY)
1700 {
1701 DEBUG1 ((" singleton Row " ID "\n", row)) ;
1702 }
1703 else
1704 {
1705 DEBUG1 ((" Row " ID " in Front " ID "\n", row, nfr)) ;
1706 }
1707 }
1708 #endif
1709
1710 /* ---------------------------------------------------------------------- */
1711 /* copy front information into Symbolic object */
1712 /* ---------------------------------------------------------------------- */
1713
1714 k = n1 ;
1715 for (i = 0 ; i < nfr ; i++)
1716 {
1717 fpivcol = Fr_npivcol [i] ;
1718 DEBUG1 (("Front " ID " k " ID " npivcol " ID " nrows " ID " ncols " ID "\n",
1719 i, k, fpivcol, Fr_nrows [i], Fr_ncols [i])) ;
1720 k += fpivcol ;
1721 /* copy Front info into Symbolic object from SW */
1722 Front_npivcol [i] = fpivcol ;
1723 Front_parent [i] = Fr_parent [i] ;
1724 }
1725
1726 /* assign empty columns to dummy placehold front nfr */
1727 DEBUG1 (("Dummy Cols in Front " ID " : " ID "\n", nfr, n_col-k)) ;
1728 Front_npivcol [nfr] = n_col - k ;
1729 Front_parent [nfr] = EMPTY ;
1730
1731 /* ---------------------------------------------------------------------- */
1732 /* find initial row permutation */
1733 /* ---------------------------------------------------------------------- */
1734
1735 /* order the singleton pivot rows */
1736 for (k = 0 ; k < n1 ; k++)
1737 {
1738 Rperm_init [k] = Rperm1 [k] ;
1739 }
1740
1741 /* determine the first row in each front (in the new row ordering) */
1742 for (i = 0 ; i < nfr ; i++)
1743 {
1744 f1rows = Front_1strow [i] ;
1745 DEBUG1 (("Front " ID " : npivcol " ID " parent " ID,
1746 i, Front_npivcol [i], Front_parent [i])) ;
1747 DEBUG1 ((" 1st rows in Front " ID " : " ID "\n", i, f1rows)) ;
1748 Front_1strow [i] = k ;
1749 k += f1rows ;
1750 }
1751
1752 /* assign empty rows to dummy placehold front nfr */
1753 DEBUG1 (("Rows in Front " ID " (dummy): " ID "\n", nfr, n_row-k)) ;
1754 Front_1strow [nfr] = k ;
1755 DEBUG1 (("nfr " ID " 1strow[nfr] " ID " nrow " ID "\n", nfr, k, n_row)) ;
1756
1757 /* Use Ci as temporary workspace for F1 */
1758 F1 = Ci ; /* [ of size nfr+1 */
1759 ASSERT (Clen >= 2*n_row + nfr+1) ;
1760
1761 for (i = 0 ; i <= nfr ; i++)
1762 {
1763 F1 [i] = Front_1strow [i] ;
1764 }
1765
1766 for (row = 0 ; row < n_row ; row++)
1767 {
1768 i = InFront [row] ;
1769 if (i != EMPTY)
1770 {
1771 newrow = F1 [i]++ ;
1772 ASSERT (newrow >= n1) ;
1773 Rperm_init [newrow] = row ;
1774 }
1775 }
1776 Rperm_init [n_row] = EMPTY ; /* unused */
1777
1778 #ifndef NDEBUG
1779 for (k = 0 ; k < n_row ; k++)
1780 {
1781 DEBUG2 (("Rperm_init [" ID "] = " ID "\n", k, Rperm_init [k])) ;
1782 }
1783 #endif
1784
1785 /* ] done using F1 */
1786
1787 /* ---------------------------------------------------------------------- */
1788 /* find the diagonal map */
1789 /* ---------------------------------------------------------------------- */
1790
1791 /* Rperm_init [newrow] = row gives the row permutation that is implied
1792 * by the column permutation, where "row" is a row index of the original
1793 * matrix A. It is not dependent on the Rperm_2by2 permutation, which
1794 * only redefines the "diagonal". Both are used to construct the
1795 * Diagonal_map. Diagonal_map only needs to be defined for
1796 * k = n1 to nn - nempty, but go ahead and define it for all of
1797 * k = 0 to nn */
1798
1799 if (prefer_diagonal)
1800 {
1801 Int *Diagonal_map ;
1802 ASSERT (n_row == n_col && nn == n_row) ;
1803 ASSERT (nempty_row == nempty_col && nempty == nempty_row) ;
1804
1805 /* allocate the Diagonal_map */
1806 Symbolic->Diagonal_map = (Int *) UMF_malloc (n_col+1, sizeof (Int)) ;
1807 Diagonal_map = Symbolic->Diagonal_map ;
1808 if (Diagonal_map == (Int *) NULL)
1809 {
1810 /* :: out of memory (diagonal map) :: */
1811 DEBUGm4 (("out of memory: Diagonal map\n")) ;
1812 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
1813 error (&Symbolic, SW) ;
1814 return (UMFPACK_ERROR_out_of_memory) ;
1815 }
1816
1817 /* use Ci as workspace to compute the inverse of Rperm_init [ */
1818 for (newrow = 0 ; newrow < nn ; newrow++)
1819 {
1820 oldrow = Rperm_init [newrow] ;
1821 ASSERT (oldrow >= 0 && oldrow < nn) ;
1822 Ci [oldrow] = newrow ;
1823 }
1824 if (strategy == UMFPACK_STRATEGY_2BY2)
1825 {
1826 ASSERT (Rperm_2by2 != (Int *) NULL) ;
1827 for (newcol = 0 ; newcol < nn ; newcol++)
1828 {
1829 oldcol = Cperm_init [newcol] ;
1830 /* 2-by-2 pivoting done in S */
1831 oldrow = Rperm_2by2 [oldcol] ;
1832 newrow = Ci [oldrow] ;
1833 Diagonal_map [newcol] = newrow ;
1834 }
1835 }
1836 else
1837 {
1838 for (newcol = 0 ; newcol < nn ; newcol++)
1839 {
1840 oldcol = Cperm_init [newcol] ;
1841 /* no 2-by-2 pivoting in S */
1842 oldrow = oldcol ;
1843 newrow = Ci [oldrow] ;
1844 Diagonal_map [newcol] = newrow ;
1845 }
1846 }
1847
1848 #ifndef NDEBUG
1849 DEBUG1 (("\nDiagonal map:\n")) ;
1850 for (newcol = 0 ; newcol < nn ; newcol++)
1851 {
1852 oldcol = Cperm_init [newcol] ;
1853 DEBUG3 (("oldcol " ID " newcol " ID ":\n", oldcol, newcol)) ;
1854 for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
1855 {
1856 Entry aij ;
1857 CLEAR (aij) ;
1858 oldrow = Ai [p] ;
1859 newrow = Ci [oldrow] ;
1860 if (Ax != (double *) NULL)
1861 {
1862 ASSIGN (aij, Ax, Az, p, SPLIT (Az)) ;
1863 }
1864 if (oldrow == oldcol)
1865 {
1866 DEBUG2 ((" old diagonal : oldcol " ID " oldrow " ID " ",
1867 oldcol, oldrow)) ;
1868 EDEBUG2 (aij) ;
1869 DEBUG2 (("\n")) ;
1870 }
1871 if (newrow == Diagonal_map [newcol])
1872 {
1873 DEBUG2 ((" MAP diagonal : newcol " ID " MAProw " ID " ",
1874 newcol, Diagonal_map [newrow])) ;
1875 EDEBUG2 (aij) ;
1876 DEBUG2 (("\n")) ;
1877 }
1878 }
1879 }
1880 #endif
1881 /* done using Ci as workspace ] */
1882
1883 }
1884
1885 /* ---------------------------------------------------------------------- */
1886 /* find the leftmost descendant of each front */
1887 /* ---------------------------------------------------------------------- */
1888
1889 for (i = 0 ; i <= nfr ; i++)
1890 {
1891 Front_leftmostdesc [i] = EMPTY ;
1892 }
1893
1894 for (i = 0 ; i < nfr ; i++)
1895 {
1896 /* start at i and walk up the tree */
1897 DEBUG2 (("Walk up front tree from " ID "\n", i)) ;
1898 j = i ;
1899 while (j != EMPTY && Front_leftmostdesc [j] == EMPTY)
1900 {
1901 DEBUG3 ((" Leftmost desc of " ID " is " ID "\n", j, i)) ;
1902 Front_leftmostdesc [j] = i ;
1903 j = Front_parent [j] ;
1904 DEBUG3 ((" go to j = " ID "\n", j)) ;
1905 }
1906 }
1907
1908 /* ---------------------------------------------------------------------- */
1909 /* find the frontal matrix chains and max frontal matrix sizes */
1910 /* ---------------------------------------------------------------------- */
1911
1912 maxnrows = 1 ; /* max # rows in any front */
1913 maxncols = 1 ; /* max # cols in any front */
1914 dmaxfrsize = 1 ; /* max frontal matrix size */
1915
1916 /* start the first chain */
1917 nchains = 0 ; /* number of chains */
1918 Chain_start [0] = 0 ; /* front 0 starts a new chain */
1919 maxrows = 1 ; /* max # rows for any front in current chain */
1920 maxcols = 1 ; /* max # cols for any front in current chain */
1921 DEBUG1 (("Constructing chains:\n")) ;
1922
1923 for (i = 0 ; i < nfr ; i++)
1924 {
1925 /* get frontal matrix info */
1926 fpivcol = Front_npivcol [i] ; /* # candidate pivot columns */
1927 fallrows = Fr_nrows [i] ; /* all rows (not just Schur comp) */
1928 fallcols = Fr_ncols [i] ; /* all cols (not just Schur comp) */
1929 parent = Front_parent [i] ; /* parent in column etree */
1930 fpiv = MIN (fpivcol, fallrows) ; /* # pivot rows and cols */
1931 maxrows = MAX (maxrows, fallrows) ;
1932 maxcols = MAX (maxcols, fallcols) ;
1933
1934 DEBUG1 (("Front: " ID ", pivcol " ID ", " ID "-by-" ID " parent " ID
1935 ", npiv " ID " Chain: maxrows " ID " maxcols " ID "\n", i, fpivcol,
1936 fallrows, fallcols, parent, fpiv, maxrows, maxcols)) ;
1937
1938 if (parent != i+1)
1939 {
1940 /* this is the end of a chain */
1941 double s ;
1942 DEBUG1 (("\nEnd of chain " ID "\n", nchains)) ;
1943
1944 /* make sure maxrows is an odd number */
1945 ASSERT (maxrows >= 0) ;
1946 if (maxrows % 2 == 0) maxrows++ ;
1947
1948 DEBUG1 (("Chain maxrows " ID " maxcols " ID "\n", maxrows, maxcols)) ;
1949
1950 Chain_maxrows [nchains] = maxrows ;
1951 Chain_maxcols [nchains] = maxcols ;
1952
1953 /* keep track of the maximum front size for all chains */
1954
1955 /* for Info only: */
1956 s = (double) maxrows * (double) maxcols ;
1957 dmaxfrsize = MAX (dmaxfrsize, s) ;
1958
1959 /* for the subsequent numerical factorization */
1960 maxnrows = MAX (maxnrows, maxrows) ;
1961 maxncols = MAX (maxncols, maxcols) ;
1962
1963 DEBUG1 (("Chain dmaxfrsize %g\n\n", dmaxfrsize)) ;
1964
1965 /* start the next chain */
1966 nchains++ ;
1967 Chain_start [nchains] = i+1 ;
1968 maxrows = 1 ;
1969 maxcols = 1 ;
1970 }
1971 }
1972
1973 /* for Info only: */
1974 dmaxfrsize = ceil (dmaxfrsize) ;
1975 DEBUGm1 (("dmaxfrsize %30.20g Int_MAX " ID "\n", dmaxfrsize, Int_MAX)) ;
1976 ASSERT (Symbolic->nchains == nchains) ;
1977
1978 /* For allocating objects in umfpack_numeric (does not include all possible
1979 * pivots, particularly pivots from prior fronts in the chain. Need to add
1980 * nb to these to get the # of columns in the L block, for example. This
1981 * is the largest row dimension and largest column dimension of any frontal
1982 * matrix. maxnrows is always odd. */
1983 Symbolic->maxnrows = maxnrows ;
1984 Symbolic->maxncols = maxncols ;
1985 DEBUGm3 (("maxnrows " ID " maxncols " ID "\n", maxnrows, maxncols)) ;
1986
1987 /* ---------------------------------------------------------------------- */
1988 /* find the initial element sizes */
1989 /* ---------------------------------------------------------------------- */
1990
1991 if (max_rdeg > dense_row_threshold)
1992 {
1993 /* there are one or more dense rows in the input matrix */
1994 /* count the number of dense rows in each column */
1995 /* use Ci as workspace for inverse of Rperm_init [ */
1996 ASSERT (Esize != (Int *) NULL) ;
1997 for (newrow = 0 ; newrow < n_row ; newrow++)
1998 {
1999 oldrow = Rperm_init [newrow] ;
2000 ASSERT (oldrow >= 0 && oldrow < nn) ;
2001 Ci [oldrow] = newrow ;
2002 }
2003 for (col = n1 ; col < n_col - nempty_col ; col++)
2004 {
2005 oldcol = Cperm_init [col] ;
2006 esize = Cdeg [oldcol] ;
2007 ASSERT (esize > 0) ;
2008 for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
2009 {
2010 oldrow = Ai [p] ;
2011 newrow = Ci [oldrow] ;
2012 if (newrow >= n1 && Rdeg [oldrow] > dense_row_threshold)
2013 {
2014 esize-- ;
2015 }
2016 }
2017 ASSERT (esize >= 0) ;
2018 Esize [col - n1] = esize ;
2019 }
2020 /* done using Ci as workspace ] */
2021 }
2022
2023 /* If there are no dense rows, then Esize [col-n1] is identical to
2024 * Cdeg [col], once Cdeg is permuted below */
2025
2026 /* ---------------------------------------------------------------------- */
2027 /* permute Cdeg and Rdeg according to initial column and row permutation */
2028 /* ---------------------------------------------------------------------- */
2029
2030 /* use Ci as workspace [ */
2031 for (k = 0 ; k < n_col ; k++)
2032 {
2033 Ci [k] = Cdeg [Cperm_init [k]] ;
2034 }
2035 for (k = 0 ; k < n_col ; k++)
2036 {
2037 Cdeg [k] = Ci [k] ;
2038 }
2039 for (k = 0 ; k < n_row ; k++)
2040 {
2041 Ci [k] = Rdeg [Rperm_init [k]] ;
2042 }
2043 for (k = 0 ; k < n_row ; k++)
2044 {
2045 Rdeg [k] = Ci [k] ;
2046 }
2047 /* done using Ci as workspace ] */
2048
2049 /* ---------------------------------------------------------------------- */
2050 /* simulate UMF_kernel_init */
2051 /* ---------------------------------------------------------------------- */
2052
2053 /* count elements and tuples at tail, LU factors of singletons, and
2054 * head and tail markers */
2055
2056 dlnz = n_inner ; /* upper limit of nz in L (incl diag) */
2057 dunz = dlnz ; /* upper limit of nz in U (incl diag) */
2058
2059 /* head marker */
2060 head_usage = 1 ;
2061 dhead_usage = 1 ;
2062
2063 /* tail markers: */
2064 tail_usage = 2 ;
2065 dtail_usage = 2 ;
2066
2067 /* allocate the Rpi and Rpx workspace for UMF_kernel_init (incl. headers) */
2068 tail_usage += UNITS (Int *, n_row+1) + UNITS (Entry *, n_row+1) + 2 ;
2069 dtail_usage += DUNITS (Int *, n_row+1) + DUNITS (Entry *, n_row+1) + 2 ;
2070 DEBUG1 (("Symbolic usage after Rpi/Rpx allocation: head " ID " tail " ID "\n",
2071 head_usage, tail_usage)) ;
2072
2073 /* LU factors for singletons, at the head of memory */
2074 for (k = 0 ; k < n1 ; k++)
2075 {
2076 lnz = Cdeg [k] - 1 ;
2077 unz = Rdeg [k] - 1 ;
2078 dlnz += lnz ;
2079 dunz += unz ;
2080 DEBUG1 (("singleton k " ID " pivrow " ID " pivcol " ID " lnz " ID " unz " ID "\n",
2081 k, Rperm_init [k], Cperm_init [k], lnz, unz)) ;
2082 head_usage += UNITS (Int, lnz) + UNITS (Entry, lnz)
2083 + UNITS (Int, unz) + UNITS (Entry, unz) ;
2084 dhead_usage += DUNITS (Int, lnz) + DUNITS (Entry, lnz)
2085 + DUNITS (Int, unz) + DUNITS (Entry, unz) ;
2086 }
2087 DEBUG1 (("Symbolic init head usage: " ID " for LU singletons\n",head_usage)) ;
2088
2089 /* column elements: */
2090 for (k = n1 ; k < n_col - nempty_col; k++)
2091 {
2092 esize = Esize ? Esize [k-n1] : Cdeg [k] ;
2093 DEBUG2 ((" esize: " ID "\n", esize)) ;
2094 ASSERT (esize >= 0) ;
2095 if (esize > 0)
2096 {
2097 tail_usage += GET_ELEMENT_SIZE (esize, 1) + 1 ;
2098 dtail_usage += DGET_ELEMENT_SIZE (esize, 1) + 1 ;
2099 }
2100 }
2101
2102 /* dense row elements */
2103 if (Esize)
2104 {
2105 Int nrow_elements = 0 ;
2106 for (k = n1 ; k < n_row - nempty_row ; k++)
2107 {
2108 rdeg = Rdeg [k] ;
2109 if (rdeg > dense_row_threshold)
2110 {
2111 tail_usage += GET_ELEMENT_SIZE (1, rdeg) + 1 ;
2112 dtail_usage += GET_ELEMENT_SIZE (1, rdeg) + 1 ;
2113 nrow_elements++ ;
2114 }
2115 }
2116 Info [UMFPACK_NDENSE_ROW] = nrow_elements ;
2117 }
2118
2119 DEBUG1 (("Symbolic usage: " ID " = head " ID " + tail " ID " after els\n",
2120 head_usage + tail_usage, head_usage, tail_usage)) ;
2121
2122 /* compute the tuple lengths */
2123 if (Esize)
2124 {
2125 /* row tuples */
2126 for (row = n1 ; row < n_row ; row++)
2127 {
2128 rdeg = Rdeg [row] ;
2129 tlen = (rdeg > dense_row_threshold) ? 1 : rdeg ;
2130 tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2131 dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2132 }
2133 /* column tuples */
2134 for (col = n1 ; col < n_col - nempty_col ; col++)
2135 {
2136 /* tlen is 1 plus the number of dense rows in this column */
2137 esize = Esize [col - n1] ;
2138 tlen = (esize > 0) + (Cdeg [col] - esize) ;
2139 tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2140 dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2141 }
2142 for ( ; col < n_col ; col++)
2143 {
2144 tail_usage += 1 + UNITS (Tuple, TUPLES (0)) ;
2145 dtail_usage += 1 + DUNITS (Tuple, TUPLES (0)) ;
2146 }
2147 }
2148 else
2149 {
2150 /* row tuples */
2151 for (row = n1 ; row < n_row ; row++)
2152 {
2153 tlen = Rdeg [row] ;
2154 tail_usage += 1 + UNITS (Tuple, TUPLES (tlen)) ;
2155 dtail_usage += 1 + DUNITS (Tuple, TUPLES (tlen)) ;
2156 }
2157 /* column tuples */
2158 for (col = n1 ; col < n_col ; col++)
2159 {
2160 tail_usage += 1 + UNITS (Tuple, TUPLES (1)) ;
2161 dtail_usage += 1 + DUNITS (Tuple, TUPLES (1)) ;
2162 }
2163 }
2164
2165 Symbolic->num_mem_init_usage = head_usage + tail_usage ;
2166 DEBUG1 (("Symbolic usage: " ID " = head " ID " + tail " ID " final\n",
2167 Symbolic->num_mem_init_usage, head_usage, tail_usage)) ;
2168
2169 ASSERT (UMF_is_permutation (Rperm_init, Ci, n_row, n_row)) ;
2170
2171 /* initial head and tail usage in Numeric->Memory */
2172 dmax_usage = dhead_usage + dtail_usage ;
2173 dmax_usage = MAX (Symbolic->num_mem_init_usage, ceil (dmax_usage)) ;
2174 Info [UMFPACK_VARIABLE_INIT_ESTIMATE] = dmax_usage ;
2175
2176 /* In case Symbolic->num_mem_init_usage overflows, keep as a double, too */
2177 Symbolic->dnum_mem_init_usage = dmax_usage ;
2178
2179 /* free the Rpi and Rpx workspace */
2180 tail_usage -= UNITS (Int *, n_row+1) + UNITS (Entry *, n_row+1) ;
2181 dtail_usage -= DUNITS (Int *, n_row+1) + DUNITS (Entry *, n_row+1) ;
2182
2183 /* ---------------------------------------------------------------------- */
2184 /* simulate UMF_kernel, assuming unsymmetric pivoting */
2185 /* ---------------------------------------------------------------------- */
2186
2187 /* Use Ci as temporary workspace for link lists [ */
2188 Link = Ci ;
2189 for (i = 0 ; i < nfr ; i++)
2190 {
2191 Link [i] = EMPTY ;
2192 }
2193
2194 flops = 0 ; /* flop count upper bound */
2195
2196 for (chain = 0 ; chain < nchains ; chain++)
2197 {
2198 double fsize ;
2199 f1 = Chain_start [chain] ;
2200 f2 = Chain_start [chain+1] - 1 ;
2201
2202 /* allocate frontal matrix working array (C, L, and U) */
2203 dr = Chain_maxrows [chain] ;
2204 dc = Chain_maxcols [chain] ;
2205 fsize =
2206 nb*nb /* LU is nb-by-nb */
2207 + dr*nb /* L is dr-by-nb */
2208 + nb*dc /* U is nb-by-dc, stored by rows */
2209 + dr*dc ; /* C is dr by dc */
2210 dtail_usage += DUNITS (Entry, fsize) ;
2211 dmax_usage = MAX (dmax_usage, dhead_usage + dtail_usage) ;
2212
2213 for (i = f1 ; i <= f2 ; i++)
2214 {
2215
2216 /* get frontal matrix info */
2217 fpivcol = Front_npivcol [i] ; /* # candidate pivot columns */
2218 fallrows = Fr_nrows [i] ; /* all rows (not just Schur comp*/
2219 fallcols = Fr_ncols [i] ; /* all cols (not just Schur comp*/
2220 parent = Front_parent [i] ; /* parent in column etree */
2221 fpiv = MIN (fpivcol, fallrows) ; /* # pivot rows and cols */
2222 f = (double) fpiv ;
2223 r = fallrows - fpiv ; /* # rows in Schur comp. */
2224 c = fallcols - fpiv ; /* # cols in Schur comp. */
2225
2226 /* assemble all children of front i in column etree */
2227 for (child = Link [i] ; child != EMPTY ; child = Link [child])
2228 {
2229 ASSERT (child >= 0 && child < i) ;
2230 ASSERT (Front_parent [child] == i) ;
2231 /* free the child element and remove it from tuple lists */
2232 cp = MIN (Front_npivcol [child], Fr_nrows [child]) ;
2233 cr = Fr_nrows [child] - cp ;
2234 cc = Fr_ncols [child] - cp ;
2235 ASSERT (cp >= 0 && cr >= 0 && cc >= 0) ;
2236 dtail_usage -= ELEMENT_SIZE (cr, cc) ;
2237
2238 }
2239
2240 /* The flop count computed here is "canonical". */
2241
2242 /* factorize the frontal matrix */
2243 flops += DIV_FLOPS * (f*r + (f-1)*f/2) /* scale pivot columns */
2244 /* f outer products: */
2245 + MULTSUB_FLOPS * (f*r*c + (r+c)*(f-1)*f/2 + (f-1)*f*(2*f-1)/6);
2246
2247 /* count nonzeros and memory usage in double precision */
2248 dlf = (f*f-f)/2 + f*r ; /* nz in L below diagonal */
2249 duf = (f*f-f)/2 + f*c ; /* nz in U above diagonal */
2250 dlnz += dlf ;
2251 dunz += duf ;
2252
2253 /* store f columns of L and f rows of U */
2254 dhead_usage +=
2255 DUNITS (Entry, dlf + duf) /* numerical values (excl diag) */
2256 + DUNITS (Int, r + c + f) ; /* indices (compressed) */
2257
2258 if (parent != EMPTY)
2259 {
2260 /* create new element and place in tuple lists */
2261 dtail_usage += ELEMENT_SIZE (r, c) ;
2262
2263 /* place in link list of parent */
2264 Link [i] = Link [parent] ;
2265 Link [parent] = i ;
2266 }
2267
2268 /* keep track of peak Numeric->Memory usage */
2269 dmax_usage = MAX (dmax_usage, dhead_usage + dtail_usage) ;
2270
2271 }
2272
2273 /* free the current frontal matrix */
2274 dtail_usage -= DUNITS (Entry, fsize) ;
2275 }
2276
2277 dhead_usage = ceil (dhead_usage) ;
2278 dmax_usage = ceil (dmax_usage) ;
2279 Symbolic->num_mem_size_est = dhead_usage ;
2280 Symbolic->num_mem_usage_est = dmax_usage ;
2281 Symbolic->lunz_bound = dlnz + dunz - n_inner ;
2282
2283 /* ] done using Ci as workspace for Link array */
2284
2285 /* ---------------------------------------------------------------------- */
2286 /* estimate total memory usage in UMFPACK_numeric */
2287 /* ---------------------------------------------------------------------- */
2288
2289 UMF_set_stats (
2290 Info,
2291 Symbolic,
2292 dmax_usage, /* estimated peak size of Numeric->Memory */
2293 dhead_usage, /* estimated final size of Numeric->Memory */
2294 flops, /* estimated "true flops" */
2295 dlnz, /* estimated nz in L */
2296 dunz, /* estimated nz in U */
2297 dmaxfrsize, /* estimated largest front size */
2298 (double) n_col, /* worst case Numeric->Upattern size */
2299 (double) n_inner, /* max possible pivots to be found */
2300 (double) maxnrows, /* estimated largest #rows in front */
2301 (double) maxncols, /* estimated largest #cols in front */
2302 TRUE, /* assume scaling is to be performed */
2303 prefer_diagonal,
2304 ESTIMATE) ;
2305
2306 /* ---------------------------------------------------------------------- */
2307
2308 #ifndef NDEBUG
2309 for (i = 0 ; i < nchains ; i++)
2310 {
2311 DEBUG2 (("Chain " ID " start " ID " end " ID " maxrows " ID " maxcols " ID "\n",
2312 i, Chain_start [i], Chain_start [i+1] - 1,
2313 Chain_maxrows [i], Chain_maxcols [i])) ;
2314 UMF_dump_chain (Chain_start [i], Fr_parent, Fr_npivcol, Fr_nrows,
2315 Fr_ncols, nfr) ;
2316 }
2317 fpivcol = 0 ;
2318 for (i = 0 ; i < nfr ; i++)
2319 {
2320 fpivcol = MAX (fpivcol, Front_npivcol [i]) ;
2321 }
2322 DEBUG0 (("Max pivot cols in any front: " ID "\n", fpivcol)) ;
2323 DEBUG1 (("Largest front: maxnrows " ID " maxncols " ID " dmaxfrsize %g\n",
2324 maxnrows, maxncols, dmaxfrsize)) ;
2325 #endif
2326
2327 /* ---------------------------------------------------------------------- */
2328 /* UMFPACK_symbolic was successful, return the object handle */
2329 /* ---------------------------------------------------------------------- */
2330
2331 Symbolic->valid = SYMBOLIC_VALID ;
2332 *SymbolicHandle = (void *) Symbolic ;
2333
2334 /* ---------------------------------------------------------------------- */
2335 /* free workspace */
2336 /* ---------------------------------------------------------------------- */
2337
2338 /* (6) The last of the workspace is free'd. The final Symbolic object
2339 * consists of 12 to 14 allocated objects. Its final total size is lies
2340 * roughly between 4*n and 13*n for a square matrix, which is all that is
2341 * left of the memory allocated by this routine. If an error occurs, the
2342 * entire Symbolic object is free'd when this routine returns (the error
2343 * return routine, below).
2344 */
2345
2346 free_work (SW) ;
2347
2348 DEBUG0 (("(3)Symbolic UMF_malloc_count - init_count = " ID "\n",
2349 UMF_malloc_count - init_count)) ;
2350 ASSERT (UMF_malloc_count == init_count + 12
2351 + (Symbolic->Esize != (Int *) NULL)
2352 + (Symbolic->Diagonal_map != (Int *) NULL)) ;
2353
2354 /* ---------------------------------------------------------------------- */
2355 /* get the time used by UMFPACK_*symbolic */
2356 /* ---------------------------------------------------------------------- */
2357
2358 umfpack_toc (stats) ;
2359 Info [UMFPACK_SYMBOLIC_WALLTIME] = stats [0] ;
2360 Info [UMFPACK_SYMBOLIC_TIME] = stats [1] ;
2361
2362 return (UMFPACK_OK) ;
2363 }
2364
2365
2366 /* ========================================================================== */
2367 /* === free_work ============================================================ */
2368 /* ========================================================================== */
2369
free_work(SWType * SW)2370 PRIVATE void free_work
2371 (
2372 SWType *SW
2373 )
2374 {
2375 if (SW)
2376 {
2377 SW->Rperm_2by2 = (Int *) UMF_free ((void *) SW->Rperm_2by2) ;
2378 SW->InvRperm1 = (Int *) UMF_free ((void *) SW->InvRperm1) ;
2379 SW->Rs = (double *) UMF_free ((void *) SW->Rs) ;
2380 SW->Si = (Int *) UMF_free ((void *) SW->Si) ;
2381 SW->Sp = (Int *) UMF_free ((void *) SW->Sp) ;
2382 SW->Ci = (Int *) UMF_free ((void *) SW->Ci) ;
2383 SW->Front_npivcol = (Int *) UMF_free ((void *) SW->Front_npivcol);
2384 SW->Front_nrows = (Int *) UMF_free ((void *) SW->Front_nrows) ;
2385 SW->Front_ncols = (Int *) UMF_free ((void *) SW->Front_ncols) ;
2386 SW->Front_parent = (Int *) UMF_free ((void *) SW->Front_parent) ;
2387 SW->Front_cols = (Int *) UMF_free ((void *) SW->Front_cols) ;
2388 SW->Cperm1 = (Int *) UMF_free ((void *) SW->Cperm1) ;
2389 SW->Rperm1 = (Int *) UMF_free ((void *) SW->Rperm1) ;
2390 SW->InFront = (Int *) UMF_free ((void *) SW->InFront) ;
2391 }
2392 }
2393
2394
2395 /* ========================================================================== */
2396 /* === error ================================================================ */
2397 /* ========================================================================== */
2398
2399 /* Error return from UMFPACK_symbolic. Free all allocated memory. */
2400
error(SymbolicType ** Symbolic,SWType * SW)2401 PRIVATE void error
2402 (
2403 SymbolicType **Symbolic,
2404 SWType *SW
2405 )
2406 {
2407
2408 free_work (SW) ;
2409 UMFPACK_free_symbolic ((void **) Symbolic) ;
2410 ASSERT (UMF_malloc_count == init_count) ;
2411 }
2412