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