|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === umf_internal.h ======================================================= */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* -------------------------------------------------------------------------- */ |
|
|
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */ |
|
|
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */ |
|
|
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */ |
|
|
9 /* -------------------------------------------------------------------------- */ |
|
|
10 |
|
|
11 /* |
|
|
12 This file is for internal use in UMFPACK itself, and should not be included |
|
|
13 in user code. Use umfpack.h instead. User-accessible file names and |
|
|
14 routine names all start with the letters "umfpack_". Non-user-accessible |
|
|
15 file names and routine names all start with "umf_". |
|
|
16 */ |
|
|
17 |
|
|
18 #ifndef _UMF_INTERNAL |
|
|
19 #define _UMF_INTERNAL |
|
|
20 |
|
|
21 /* -------------------------------------------------------------------------- */ |
|
|
22 /* ANSI standard include files */ |
|
|
23 /* -------------------------------------------------------------------------- */ |
|
|
24 |
|
|
25 /* from float.h: DBL_EPSILON */ |
|
|
26 #include <float.h> |
|
|
27 |
|
|
28 /* from string.h: strcmp */ |
|
|
29 #include <string.h> |
|
|
30 |
|
|
31 /* when debugging, assert.h and the assert macro are used (see umf_dump.h) */ |
|
|
32 |
|
|
33 /* -------------------------------------------------------------------------- */ |
|
|
34 /* Architecture */ |
|
|
35 /* -------------------------------------------------------------------------- */ |
|
|
36 |
|
|
37 #if defined (__sun) || defined (MSOL2) || defined (ARCH_SOL2) |
|
|
38 #define UMF_SOL2 |
|
|
39 #define UMFPACK_ARCHITECTURE "Sun Solaris" |
|
|
40 |
|
|
41 #elif defined (__sgi) || defined (MSGI) || defined (ARCH_SGI) |
|
|
42 #define UMF_SGI |
|
|
43 #define UMFPACK_ARCHITECTURE "SGI Irix" |
|
|
44 |
|
|
45 #elif defined (__linux) || defined (MGLNX86) || defined (ARCH_GLNX86) |
|
|
46 #define UMF_LINUX |
|
|
47 #define UMFPACK_ARCHITECTURE "Linux" |
|
|
48 |
|
|
49 #elif defined (_AIX) || defined (MIBM_RS) || defined (ARCH_IBM_RS) |
|
|
50 #define UMF_AIX |
|
|
51 #define UMFPACK_ARCHITECTURE "IBM AIX" |
|
|
52 |
|
|
53 #elif defined (__alpha) || defined (MALPHA) || defined (ARCH_ALPHA) |
|
|
54 #define UMF_ALPHA |
|
|
55 #define UMFPACK_ARCHITECTURE "Compaq Alpha" |
|
|
56 |
|
|
57 #elif defined (__WIN32) || defined (_WIN32) || defined (_win32) || defined (__win32) || defined (WIN32) |
|
|
58 #define UMF_WINDOWS |
|
|
59 #define UMFPACK_ARCHITECTURE "Microsoft Windows" |
|
|
60 |
|
|
61 #elif defined (__hppa) || defined (__hpux) || defined (MHPUX) || defined (ARCH_HPUX) |
|
|
62 #define UMF_HP |
|
|
63 #define UMFPACK_ARCHITECTURE "HP Unix" |
|
|
64 |
|
|
65 #elif defined (__hp700) || defined (MHP700) || defined (ARCH_HP700) |
|
|
66 #define UMF_HP |
|
|
67 #define UMFPACK_ARCHITECTURE "HP 700 Unix" |
|
|
68 |
|
|
69 #else |
|
|
70 /* If the architecture is unknown, and you call the BLAS, you may need to */ |
|
|
71 /* define BLAS_BY_VALUE, BLAS_NO_UNDERSCORE, and/or BLAS_CHAR_ARG yourself. */ |
|
|
72 #define UMFPACK_ARCHITECTURE "unknown" |
|
|
73 #endif |
|
|
74 |
|
|
75 |
|
|
76 /* -------------------------------------------------------------------------- */ |
|
|
77 /* basic definitions (see also amd_internal.h) */ |
|
|
78 /* -------------------------------------------------------------------------- */ |
|
|
79 |
|
|
80 #define ONES_COMPLEMENT(r) (-(r)-1) |
|
|
81 |
|
|
82 /* -------------------------------------------------------------------------- */ |
|
|
83 /* AMD include file */ |
|
|
84 /* -------------------------------------------------------------------------- */ |
|
|
85 |
|
|
86 /* stdio.h, stdlib.h, limits.h, and math.h, NDEBUG definition, |
|
|
87 * assert.h, and MATLAB include files */ |
|
|
88 #include "amd_internal.h" |
|
|
89 |
|
|
90 /* -------------------------------------------------------------------------- */ |
|
|
91 /* Real/complex and int/long definitions, double relops */ |
|
|
92 /* -------------------------------------------------------------------------- */ |
|
|
93 |
|
|
94 #include "umf_version.h" |
|
|
95 |
|
|
96 /* -------------------------------------------------------------------------- */ |
|
|
97 /* Compile-time configurations */ |
|
|
98 /* -------------------------------------------------------------------------- */ |
|
|
99 |
|
|
100 #include "umf_config.h" |
|
|
101 |
|
|
102 /* -------------------------------------------------------------------------- */ |
|
|
103 /* umfpack include file */ |
|
|
104 /* -------------------------------------------------------------------------- */ |
|
|
105 |
|
|
106 #include "umfpack.h" |
|
|
107 |
|
|
108 /* -------------------------------------------------------------------------- */ |
|
|
109 /* for contents of Info. This must correlate with umfpack.h */ |
|
|
110 /* -------------------------------------------------------------------------- */ |
|
|
111 |
|
|
112 #define ESTIMATE (UMFPACK_NUMERIC_SIZE_ESTIMATE - UMFPACK_NUMERIC_SIZE) |
|
|
113 #define ACTUAL 0 |
|
|
114 |
|
|
115 /* -------------------------------------------------------------------------- */ |
|
|
116 /* get a parameter from the Control array */ |
|
|
117 /* -------------------------------------------------------------------------- */ |
|
|
118 |
|
|
119 #define GET_CONTROL(i,default) \ |
|
|
120 ((Control != (double *) NULL) ? \ |
|
|
121 (SCALAR_IS_NAN (Control [i]) ? default : Control [i]) \ |
|
|
122 : default) |
|
|
123 |
|
|
124 /* -------------------------------------------------------------------------- */ |
|
|
125 /* for clearing the external degree counters */ |
|
|
126 /* -------------------------------------------------------------------------- */ |
|
|
127 |
|
|
128 #define MAX_MARK(n) Int_MAX - (2*(n)+1) |
|
|
129 |
|
|
130 /* -------------------------------------------------------------------------- */ |
|
|
131 /* convert number of Units to MBytes */ |
|
|
132 /* -------------------------------------------------------------------------- */ |
|
|
133 |
|
|
134 #define MBYTES(units) (((units) * sizeof (Unit)) / 1048576.0) |
|
|
135 |
|
|
136 /* -------------------------------------------------------------------------- */ |
|
|
137 /* dense row/column macro */ |
|
|
138 /* -------------------------------------------------------------------------- */ |
|
|
139 |
|
|
140 /* In order for a row or column to be treated as "dense", it must have more */ |
|
|
141 /* entries than the value returned by this macro. n is the dimension of the */ |
|
|
142 /* matrix, and alpha is the dense row/column control parameter. */ |
|
|
143 |
|
|
144 /* Note: this is not defined if alpha is NaN or Inf: */ |
|
|
145 #define UMFPACK_DENSE_DEGREE_THRESHOLD(alpha,n) \ |
|
|
146 ((Int) MAX (16.0, (alpha) * 16.0 * sqrt ((double) (n)))) |
|
|
147 |
|
|
148 /* -------------------------------------------------------------------------- */ |
|
|
149 /* PRINTF */ |
|
|
150 /* -------------------------------------------------------------------------- */ |
|
|
151 |
|
|
152 #define PRINTFk(k,params) { if (prl >= (k)) { PRINTF (params) ; } } |
|
|
153 #define PRINTF1(params) PRINTFk (1, params) |
|
|
154 #define PRINTF2(params) PRINTFk (2, params) |
|
|
155 #define PRINTF3(params) PRINTFk (3, params) |
|
|
156 #define PRINTF4(params) PRINTFk (4, params) |
|
|
157 #define PRINTF5(params) PRINTFk (5, params) |
|
|
158 #define PRINTF6(params) PRINTFk (6, params) |
|
|
159 |
|
|
160 /* -------------------------------------------------------------------------- */ |
|
|
161 /* Fixed control parameters */ |
|
|
162 /* -------------------------------------------------------------------------- */ |
|
|
163 |
|
|
164 /* maximum number of columns to consider at one time, in a single front */ |
|
|
165 #define MAX_CANDIDATES 128 |
|
|
166 |
|
|
167 /* reduce Numeric->Memory request by this ratio, if allocation fails */ |
|
|
168 #define UMF_REALLOC_REDUCTION (0.95) |
|
|
169 |
|
|
170 /* increase Numeric->Memory request by this ratio, if we need more */ |
|
|
171 #define UMF_REALLOC_INCREASE (1.2) |
|
|
172 |
|
|
173 /* increase the dimensions of the current frontal matrix by this factor |
|
|
174 * when it needs to grow. */ |
|
|
175 #define UMF_FRONTAL_GROWTH (1.2) |
|
|
176 |
|
|
177 /* largest BLAS block size permitted */ |
|
|
178 #define MAXNB 64 |
|
|
179 |
|
|
180 /* if abs (y) < RECIPROCAL_TOLERANCE, then compute x/y. Otherwise x*(1/y). |
|
|
181 * Ignored if NRECIPROCAL is defined */ |
|
|
182 #define RECIPROCAL_TOLERANCE 1e-12 |
|
|
183 |
|
|
184 /* -------------------------------------------------------------------------- */ |
|
|
185 /* Memory allocator */ |
|
|
186 /* -------------------------------------------------------------------------- */ |
|
|
187 |
|
|
188 /* The MATLAB mexFunction uses MATLAB's memory manager, while the C-callable |
|
|
189 * AMD library uses the ANSI C malloc, free, and realloc routines. To use |
|
|
190 * the mx* memory allocation routines, use -DNUTIL when compiling. |
|
|
191 */ |
|
|
192 |
|
|
193 #undef ALLOCATE |
|
|
194 #undef FREE |
|
|
195 #undef REALLOC |
|
|
196 |
|
|
197 #ifdef MATLAB_MEX_FILE |
|
|
198 |
|
|
199 #ifdef NUTIL |
|
|
200 |
|
|
201 /* These functions simply terminate the mexFunction if they fail to allocate |
|
|
202 * memory. That's too restrictive for UMFPACK. */ |
|
|
203 #define ALLOCATE mxMalloc |
|
|
204 #define FREE mxFree |
|
|
205 #define REALLOCATE mxRealloc |
|
|
206 |
|
|
207 #else |
|
|
208 |
|
|
209 /* Use internal MATLAB memory allocation routines, used by built-in MATLAB |
|
|
210 * functions. These are not documented, but are available for use. Their |
|
|
211 * prototypes are in util.h, but that file is not provided to the MATLAB user. |
|
|
212 * The advantage of using these routines is that they return NULL if out of |
|
|
213 * memory, instead of terminating the mexFunction. UMFPACK attempts to allocate |
|
|
214 * extra space for "elbow room", and then reduces its request if the memory is |
|
|
215 * not available. That strategy doesn't work with the mx* routines. |
|
|
216 */ |
|
|
217 void *utMalloc (size_t size) ; |
|
|
218 void utFree (void *p) ; |
|
|
219 void *utRealloc (void *p, size_t size) ; |
|
|
220 #define ALLOCATE utMalloc |
|
|
221 #define FREE utFree |
|
|
222 #define REALLOCATE utRealloc |
|
|
223 |
|
|
224 #endif |
|
|
225 #else |
|
|
226 #ifdef MATHWORKS |
|
|
227 |
|
|
228 /* Compiling as a built-in routine. Since out-of-memory conditions are checked |
|
|
229 * after every allocation, we can use ut* routines here. */ |
|
|
230 #define ALLOCATE utMalloc |
|
|
231 #define FREE utFree |
|
|
232 #define REALLOCATE utRealloc |
|
|
233 |
|
|
234 #else |
|
|
235 |
|
|
236 /* use the ANSI C memory allocation routines */ |
|
|
237 #define ALLOCATE malloc |
|
|
238 #define FREE free |
|
|
239 #define REALLOCATE realloc |
|
|
240 |
|
|
241 #endif |
|
|
242 #endif |
|
|
243 |
|
|
244 /* -------------------------------------------------------------------------- */ |
|
|
245 /* Memory space definitions */ |
|
|
246 /* -------------------------------------------------------------------------- */ |
|
|
247 |
|
|
248 /* for memory alignment - assume double has worst case alignment */ |
|
|
249 typedef double Align ; |
|
|
250 |
|
|
251 /* get number of bytes required to hold n items of a type: */ |
|
|
252 /* note that this will not overflow, because sizeof (type) is always */ |
|
|
253 /* greater than or equal to sizeof (Int) >= 2 */ |
|
|
254 #define BYTES(type,n) (sizeof (type) * (n)) |
|
|
255 |
|
|
256 /* ceiling of (b/u). Assumes b >= 0 and u > 0 */ |
|
|
257 #define CEILING(b,u) (((b) + (u) - 1) / (u)) |
|
|
258 |
|
|
259 /* get number of Units required to hold n items of a type: */ |
|
|
260 #define UNITS(type,n) (CEILING (BYTES (type, n), sizeof (Unit))) |
|
|
261 |
|
|
262 /* same as DUNITS, but use double instead of int to avoid overflow */ |
|
|
263 #define DUNITS(type,n) (ceil (BYTES (type, (double) n) / sizeof (Unit))) |
|
|
264 |
|
|
265 union Unit_union |
|
|
266 { /* memory is allocated in multiples of Unit */ |
|
|
267 struct |
|
|
268 { |
|
|
269 Int |
|
|
270 size, /* size, in Units, of the block, excl. header block */ |
|
|
271 /* size >= 0: block is in use */ |
|
|
272 /* size < 0: block is free, of |size| Units */ |
|
|
273 prevsize ; /* size, in Units, of preceding block in S->Memory */ |
|
|
274 /* during garbage_collection, prevsize is set to -e-1 */ |
|
|
275 /* for element e, or positive (and thus a free block) */ |
|
|
276 /* otherwise */ |
|
|
277 } header ; /* block header */ |
|
|
278 Align xxxxxx ; /* force alignment of blocks (xxxxxx is never used) */ |
|
|
279 } ; |
|
|
280 |
|
|
281 typedef union Unit_union Unit ; |
|
|
282 |
|
|
283 /* get the size of an allocated block */ |
|
|
284 #define GET_BLOCK_SIZE(p) (((p)-1)->header.size) |
|
|
285 |
|
|
286 /* -------------------------------------------------------------------------- */ |
|
|
287 /* Numeric */ |
|
|
288 /* -------------------------------------------------------------------------- */ |
|
|
289 |
|
|
290 /* |
|
|
291 NUMERIC_VALID and SYMBOLIC_VALID: |
|
|
292 The different values of SYBOLIC_VALID and NUMERIC_VALID are chosen as a |
|
|
293 first defense against corrupted *Symbolic or *Numeric pointers passed to an |
|
|
294 UMFPACK routine. They also ensure that the objects are used only by the |
|
|
295 same version that created them (umfpack_di_*, umfpack_dl_*, umfpack_zi_*, |
|
|
296 or umfpack_zl_*). The values have also been changed since prior releases of |
|
|
297 the code to ensure that all routines that operate on the objects are of the |
|
|
298 same release. The values themselves are purely arbitrary. The are less |
|
|
299 than the ANSI C required minimums of INT_MAX and LONG_MAX, respectively. |
|
|
300 */ |
|
|
301 |
|
|
302 #ifdef DINT |
|
|
303 #define NUMERIC_VALID 15977 |
|
|
304 #define SYMBOLIC_VALID 41937 |
|
|
305 #endif |
|
|
306 #ifdef DLONG |
|
|
307 #define NUMERIC_VALID 399789720 |
|
|
308 #define SYMBOLIC_VALID 399192713 |
|
|
309 #endif |
|
|
310 #ifdef ZINT |
|
|
311 #define NUMERIC_VALID 17957 |
|
|
312 #define SYMBOLIC_VALID 40927 |
|
|
313 #endif |
|
|
314 #ifdef ZLONG |
|
|
315 #define NUMERIC_VALID 129987754 |
|
|
316 #define SYMBOLIC_VALID 110291734 |
|
|
317 #endif |
|
|
318 |
|
|
319 typedef struct /* NumericType */ |
|
|
320 { |
|
|
321 double |
|
|
322 flops, /* "true" flop count */ |
|
|
323 relpt, /* relative pivot tolerance used */ |
|
|
324 relpt2, /* relative pivot tolerance used for sym. */ |
|
|
325 droptol, |
|
|
326 alloc_init, /* initial allocation of Numeric->memory */ |
|
|
327 front_alloc_init, /* frontal matrix allocation parameter */ |
|
|
328 rsmin, /* smallest row sum */ |
|
|
329 rsmax, /* largest row sum */ |
|
|
330 min_udiag, /* smallest abs value on diagonal of D */ |
|
|
331 max_udiag, /* smallest abs value on diagonal of D */ |
|
|
332 rcond ; /* min (D) / max (D) */ |
|
|
333 |
|
|
334 Int |
|
|
335 scale ; |
|
|
336 |
|
|
337 Int valid ; /* set to NUMERIC_VALID, for validity check */ |
|
|
338 |
|
|
339 /* Memory space for A and LU factors */ |
|
|
340 Unit |
|
|
341 *Memory ; /* working memory for A and LU factors */ |
|
|
342 Int |
|
|
343 ihead, /* pointer to tail of LU factors, in Numeric->Memory */ |
|
|
344 itail, /* pointer to top of elements & tuples, */ |
|
|
345 /* in Numeric->Memory */ |
|
|
346 ibig, /* pointer to largest free block seen in tail */ |
|
|
347 size ; /* size of Memory, in Units */ |
|
|
348 |
|
|
349 Int |
|
|
350 *Rperm, /* pointer to row perm array, size: n+1 */ |
|
|
351 /* after UMF_kernel: Rperm [new] = old */ |
|
|
352 /* during UMF_kernel: Rperm [old] = new */ |
|
|
353 *Cperm, /* pointer to col perm array, size: n+1 */ |
|
|
354 /* after UMF_kernel: Cperm [new] = old */ |
|
|
355 /* during UMF_kernel: Cperm [old] = new */ |
|
|
356 |
|
|
357 *Upos, /* see UMFPACK_get_numeric for a description */ |
|
|
358 *Lpos, |
|
|
359 *Lip, |
|
|
360 *Lilen, |
|
|
361 *Uip, |
|
|
362 *Uilen, |
|
|
363 *Upattern ; /* pattern of last row of U (if singular) */ |
|
|
364 |
|
|
365 Int |
|
|
366 ulen, /* length of Upattern */ |
|
|
367 npiv, /* number of structural pivots found (sprank approx) */ |
|
|
368 nnzpiv ; /* number of numerical (nonzero) pivots found */ |
|
|
369 |
|
|
370 Entry |
|
|
371 *D ; /* D [i] is the diagonal entry of U */ |
|
|
372 |
|
|
373 Int do_recip ; |
|
|
374 double *Rs ; /* scale factors for the rows of A and b */ |
|
|
375 /* do_recip FALSE: Divide row i by Rs [i] */ |
|
|
376 /* do_recip TRUE: Multiply row i by Rs [i] */ |
|
|
377 |
|
|
378 Int |
|
|
379 n_row, n_col, /* A is n_row-by-n_row */ |
|
|
380 n1 ; /* number of singletons */ |
|
|
381 |
|
|
382 /* for information only: */ |
|
|
383 Int |
|
|
384 tail_usage, /* amount of memory allocated in tail */ |
|
|
385 /* head_usage is Numeric->ihead */ |
|
|
386 init_usage, /* memory usage just after UMF_kernel_init */ |
|
|
387 max_usage, /* peak memory usage (excludes internal and external */ |
|
|
388 /* fragmentation in the tail) */ |
|
|
389 ngarbage, /* number of garbage collections performed */ |
|
|
390 nrealloc, /* number of reallocations performed */ |
|
|
391 ncostly, /* number of costly reallocations performed */ |
|
|
392 isize, /* size of integer pattern of L and U */ |
|
|
393 nLentries, /* number of entries in L, excluding diagonal */ |
|
|
394 nUentries, /* number of entries in U, including diagonal */ |
|
|
395 /* Some entries may be numerically zero. */ |
|
|
396 lnz, /* number of nonzero entries in L, excl. diagonal */ |
|
|
397 all_lnz, /* lnz plus entries dropped from L */ |
|
|
398 unz, /* number of nonzero entries in U, excl. diagonal */ |
|
|
399 all_unz, /* unz plus entries dropped form U */ |
|
|
400 maxfrsize ; /* largest actual front size */ |
|
|
401 |
|
|
402 Int maxnrows, maxncols ; /* not the same as Symbolic->maxnrows/cols* */ |
|
|
403 |
|
|
404 } NumericType ; |
|
|
405 |
|
|
406 |
|
|
407 |
|
|
408 /* -------------------------------------------------------------------------- */ |
|
|
409 /* Element tuples for connecting elements together in a matrix */ |
|
|
410 /* -------------------------------------------------------------------------- */ |
|
|
411 |
|
|
412 typedef struct /* Tuple */ |
|
|
413 { |
|
|
414 /* The (e,f) tuples for the element lists */ |
|
|
415 Int |
|
|
416 e, /* element */ |
|
|
417 f ; /* contribution to the row/col appears at this offset */ |
|
|
418 |
|
|
419 } Tuple ; |
|
|
420 |
|
|
421 #define TUPLES(t) MAX (4, (t) + 1) |
|
|
422 |
|
|
423 /* Col_degree is aliased with Cperm, and Row_degree with Rperm */ |
|
|
424 #define NON_PIVOTAL_COL(col) (Col_degree [col] >= 0) |
|
|
425 #define NON_PIVOTAL_ROW(row) (Row_degree [row] >= 0) |
|
|
426 |
|
|
427 /* -------------------------------------------------------------------------- */ |
|
|
428 /* An element */ |
|
|
429 /* -------------------------------------------------------------------------- */ |
|
|
430 |
|
|
431 typedef struct /* Element */ |
|
|
432 { |
|
|
433 Int |
|
|
434 |
|
|
435 cdeg, /* external column degree + cdeg0 offset */ |
|
|
436 rdeg, /* external row degree + rdeg0 offset */ |
|
|
437 nrowsleft, /* number of rows remaining */ |
|
|
438 ncolsleft, /* number of columns remaining */ |
|
|
439 nrows, /* number of rows */ |
|
|
440 ncols, /* number of columns */ |
|
|
441 next ; /* for list link of sons, used during assembly only */ |
|
|
442 |
|
|
443 /* followed in memory by: |
|
|
444 Int |
|
|
445 col [0..ncols-1], column indices of this element |
|
|
446 row [0..nrows-1] ; row indices of this element |
|
|
447 Entry (suitably aligned, see macro below) |
|
|
448 C [0...nrows-1, 0...ncols-1] ; |
|
|
449 size of C is nrows*ncols Entry's |
|
|
450 */ |
|
|
451 |
|
|
452 } Element ; |
|
|
453 |
|
|
454 /* macros for computing pointers to row/col indices, and contribution block: */ |
|
|
455 |
|
|
456 #define GET_ELEMENT_SIZE(nr,nc) \ |
|
|
457 (UNITS (Element, 1) + UNITS (Int, (nc) + (nr)) + UNITS (Entry, (nc) * (nr))) |
|
|
458 |
|
|
459 #define DGET_ELEMENT_SIZE(nr,nc) \ |
|
|
460 (DUNITS (Element, 1) + DUNITS (Int, (nc) + (nr)) + DUNITS (Entry, (nc) * (nr))) |
|
|
461 |
|
|
462 #define GET_ELEMENT_COLS(ep,p,Cols) { \ |
|
|
463 ASSERT (p != (Unit *) NULL) ; \ |
|
|
464 ASSERT (p >= Numeric->Memory + Numeric->itail) ; \ |
|
|
465 ASSERT (p <= Numeric->Memory + Numeric->size) ; \ |
|
|
466 ep = (Element *) p ; \ |
|
|
467 p += UNITS (Element, 1) ; \ |
|
|
468 Cols = (Int *) p ; \ |
|
|
469 } |
|
|
470 |
|
|
471 #define GET_ELEMENT_PATTERN(ep,p,Cols,Rows,ncm) { \ |
|
|
472 GET_ELEMENT_COLS (ep, p, Cols) ; \ |
|
|
473 ncm = ep->ncols ; \ |
|
|
474 Rows = Cols + ncm ; \ |
|
|
475 } |
|
|
476 |
|
|
477 #define GET_ELEMENT(ep,p,Cols,Rows,ncm,nrm,C) { \ |
|
|
478 GET_ELEMENT_PATTERN (ep, p, Cols, Rows, ncm) ; \ |
|
|
479 nrm = ep->nrows ; \ |
|
|
480 p += UNITS (Int, ncm + nrm) ; \ |
|
|
481 C = (Entry *) p ; \ |
|
|
482 } |
|
|
483 |
|
|
484 /* -------------------------------------------------------------------------- */ |
|
|
485 /* Work data structure */ |
|
|
486 /* -------------------------------------------------------------------------- */ |
|
|
487 |
|
|
488 /* |
|
|
489 This data structure holds items needed only during factorization. |
|
|
490 All of this is freed when UMFPACK_numeric completes. Note that some of |
|
|
491 it is stored in the tail end of Numeric->S (namely, the Tuples and the |
|
|
492 Elements). |
|
|
493 */ |
|
|
494 |
|
|
495 typedef struct /* WorkType */ |
|
|
496 { |
|
|
497 |
|
|
498 /* ---------------------------------------------------------------------- */ |
|
|
499 /* information about each row and col of A */ |
|
|
500 /* ---------------------------------------------------------------------- */ |
|
|
501 |
|
|
502 /* |
|
|
503 Row_tuples: pointer to tuple list (alias with Numeric->Uip) |
|
|
504 Row_tlen: number of tuples (alias with Numeric->Uilen) |
|
|
505 Col_tuples: pointer to tuple list (alias with Numeric->Lip) |
|
|
506 Col_tlen: number of tuples (alias with Numeric->Lilen) |
|
|
507 Row_degree: degree of the row or column (alias Numeric->Rperm) |
|
|
508 Col_degree: degree of the row or column (alias Numeric->Cperm) |
|
|
509 |
|
|
510 The Row_degree and Col_degree are MATLAB-style colmmd approximations, |
|
|
511 are equal to the sum of the sizes of the elements (contribution blocks) |
|
|
512 in each row and column. They are maintained when elements are created |
|
|
513 and assembled. They are used only during the pivot row and column |
|
|
514 search. They are not needed to represent the pattern of the remaining |
|
|
515 matrix. |
|
|
516 */ |
|
|
517 |
|
|
518 /* ---------------------------------------------------------------------- */ |
|
|
519 /* information about each element */ |
|
|
520 /* ---------------------------------------------------------------------- */ |
|
|
521 |
|
|
522 Int *E ; /* E [0 .. Work->elen-1] element "pointers" */ |
|
|
523 /* (offsets in Numeric->Memory) */ |
|
|
524 |
|
|
525 /* ---------------------------------------------------------------------- */ |
|
|
526 /* generic workspace */ |
|
|
527 /* ---------------------------------------------------------------------- */ |
|
|
528 |
|
|
529 Entry *Wx, *Wy ; /* each of size maxnrows+1 */ |
|
|
530 |
|
|
531 Int /* Sizes: nn = MAX (n_row, n_col) */ |
|
|
532 *Wp, /* nn+1 */ |
|
|
533 *Wrp, /* n_col+1 */ |
|
|
534 *Wm, /* maxnrows+1 */ |
|
|
535 *Wio, /* maxncols+1 */ |
|
|
536 *Woi, /* maxncols+1 */ |
|
|
537 *Woo, /* MAX (maxnrows,maxncols)+1 */ |
|
|
538 *Wrow, /* pointer to Fcols, Wio, or Woi */ |
|
|
539 *NewRows, /* list of rows to scan */ |
|
|
540 *NewCols ; /* list of cols to scan */ |
|
|
541 |
|
|
542 /* ---------------------------------------------------------------------- */ |
|
|
543 |
|
|
544 Int |
|
|
545 *Lpattern, /* pattern of column of L, for one Lchain */ |
|
|
546 *Upattern, /* pattern of row of U, for one Uchain */ |
|
|
547 ulen, llen ; /* length of Upattern and Lpattern */ |
|
|
548 |
|
|
549 Int |
|
|
550 *Diagonal_map, /* used for symmetric pivoting, of size nn+1 */ |
|
|
551 *Diagonal_imap ;/* used for symmetric pivoting, of size nn+1 */ |
|
|
552 |
|
|
553 /* ---------------------------------------------------------------------- */ |
|
|
554 |
|
|
555 Int |
|
|
556 n_row, n_col, /* matrix is n_row-by-n_col */ |
|
|
557 nz, /* nonzeros in the elements for this matrix */ |
|
|
558 n1, /* number of row and col singletons */ |
|
|
559 elen, /* max possible number of elements */ |
|
|
560 npiv, /* number of pivot rows and columns so far */ |
|
|
561 ndiscard, /* number of discarded pivot columns */ |
|
|
562 Wrpflag, |
|
|
563 nel, /* elements in use are in the range 1..nel */ |
|
|
564 noff_diagonal, |
|
|
565 prior_element, |
|
|
566 rdeg0, cdeg0, |
|
|
567 rrdeg, ccdeg, |
|
|
568 Candidates [MAX_CANDIDATES], /* current candidate pivot columns */ |
|
|
569 nCandidates, /* number of candidates in Candidate set */ |
|
|
570 ksuper, |
|
|
571 firstsuper, |
|
|
572 jsuper, |
|
|
573 ncand, /* number of candidates (some not in Candidates[ ]) */ |
|
|
574 nextcand, /* next candidate to place in Candidate search set */ |
|
|
575 lo, |
|
|
576 hi, |
|
|
577 pivrow, /* current pivot row */ |
|
|
578 pivcol, /* current pivot column */ |
|
|
579 do_extend, /* true if the next pivot extends the current front */ |
|
|
580 do_update, /* true if update should be applied */ |
|
|
581 nforced, /* number of forced updates because of frontal growth */ |
|
|
582 any_skip, |
|
|
583 do_scan2row, |
|
|
584 do_scan2col, |
|
|
585 do_grow, |
|
|
586 pivot_case, |
|
|
587 frontid, /* id of current frontal matrix */ |
|
|
588 nfr ; /* number of frontal matrices */ |
|
|
589 |
|
|
590 /* ---------------------------------------------------------------------- */ |
|
|
591 /* For row-merge tree */ |
|
|
592 /* ---------------------------------------------------------------------- */ |
|
|
593 |
|
|
594 Int |
|
|
595 *Front_new1strow ; |
|
|
596 |
|
|
597 /* ---------------------------------------------------------------------- */ |
|
|
598 /* current frontal matrix, F */ |
|
|
599 /* ---------------------------------------------------------------------- */ |
|
|
600 |
|
|
601 Int Pivrow [MAXNB], |
|
|
602 Pivcol [MAXNB] ; |
|
|
603 |
|
|
604 Entry |
|
|
605 *Flublock, /* LU block, nb-by-nb */ |
|
|
606 *Flblock, /* L block, fnr_curr-by-nb */ |
|
|
607 *Fublock, /* U block, nb-by-fnc_curr, or U' fnc_curr-by-nb */ |
|
|
608 *Fcblock ; /* C block, fnr_curr-by-fnc_curr */ |
|
|
609 |
|
|
610 Int |
|
|
611 *Frows, /* Frows [0.. ]: row indices of F */ |
|
|
612 |
|
|
613 *Fcols, /* Fcols [0.. ]: column indices of F */ |
|
|
614 |
|
|
615 *Frpos, /* position of row indices in F, or -1 if not present */ |
|
|
616 /* if Frows[i] == row, then Frpos[row] == i */ |
|
|
617 |
|
|
618 *Fcpos, /* position of col indices in F, or -1 if not present */ |
|
|
619 /* if Fcols[j] == col, then */ |
|
|
620 /* Fcpos[col] == j*Work->fnr_curr */ |
|
|
621 |
|
|
622 fnrows, /* number of rows in contribution block in F */ |
|
|
623 fncols, /* number of columns in contribution block in F */ |
|
|
624 fnr_curr, /* maximum # of rows in F (leading dimension) */ |
|
|
625 fnc_curr, /* maximum # of columns in F */ |
|
|
626 fcurr_size, /* current size of F */ |
|
|
627 fnrows_max, /* max possible column-dimension (max # of rows) of F */ |
|
|
628 fncols_max, /* max possible row-dimension (max # of columns) of F */ |
|
|
629 nb, |
|
|
630 fnpiv, /* number of pivots in F */ |
|
|
631 fnzeros, /* number of explicit zero entries in LU block */ |
|
|
632 fscan_row, /* where to start scanning rows of F in UMF_assemble */ |
|
|
633 fscan_col, /* where to start scanning cols of F in UMF_assemble */ |
|
|
634 fnrows_new, /* number of new row indices in F after pivot added */ |
|
|
635 fncols_new, /* number of new col indices in F after pivot added */ |
|
|
636 pivrow_in_front, /* true if current pivot row in Frows */ |
|
|
637 pivcol_in_front ; /* true if current pivot column in Fcols */ |
|
|
638 |
|
|
639 /* ---------------------------------------------------------------------- |
|
|
640 * Current frontal matrix |
|
|
641 * ---------------------------------------------------------------------- |
|
|
642 * The current frontal matrix is held as a single block of memory allocated |
|
|
643 * from the "tail" end of Numeric->Memory. It is subdivided into four |
|
|
644 * parts: an LU block, an L block, a U block, and a C block. |
|
|
645 * |
|
|
646 * Let k = fnpiv, r = fnrows, and c = fncols for the following discussion. |
|
|
647 * Let dr = fnr_curr and dc = fnc_curr. Note that r <= dr and c <= dc. |
|
|
648 * |
|
|
649 * The LU block is of dimension nb-by-nb. The first k-by-k part holds the |
|
|
650 * "diagonal" part of the LU factors for these k pivot rows and columns. |
|
|
651 * The k pivot row and column indices in this part are Pivrow [0..k-1] and |
|
|
652 * Pivcol [0..k-1], respectively. |
|
|
653 * |
|
|
654 * The L block is of dimension dr-by-nb. It holds the k pivot columns, |
|
|
655 * except for the leading k-by-k part in the LU block. Only the leading |
|
|
656 * r-by-k part is in use. |
|
|
657 * |
|
|
658 * The U block is of dimension dc-by-nb. It holds the k pivot rows, |
|
|
659 * except for the leading k-by-k part in the LU block. It is stored in |
|
|
660 * row-oriented form. Only the leading c-by-k part is in use. |
|
|
661 * |
|
|
662 * The C block is of dimension dr-by-dc. It holds the current contribution |
|
|
663 * block. Only the leading r-by-c part is in use. The column indices in |
|
|
664 * the C block are Fcols [0..c-1], and the row indices are Frows [0..r-1]. |
|
|
665 * |
|
|
666 * dr is always odd, to avoid bad cache behavior. |
|
|
667 */ |
|
|
668 |
|
|
669 } WorkType ; |
|
|
670 |
|
|
671 |
|
|
672 /* -------------------------------------------------------------------------- */ |
|
|
673 /* Symbolic */ |
|
|
674 /* -------------------------------------------------------------------------- */ |
|
|
675 |
|
|
676 /* |
|
|
677 This is is constructed by UMFPACK_symbolic, and is needed by UMFPACK_numeric |
|
|
678 to factor the matrix. |
|
|
679 */ |
|
|
680 |
|
|
681 typedef struct /* SymbolicType */ |
|
|
682 { |
|
|
683 |
|
|
684 double |
|
|
685 num_mem_usage_est, /* estimated max Numeric->Memory size */ |
|
|
686 num_mem_size_est, /* estimated final Numeric->Memory size */ |
|
|
687 peak_sym_usage, /* peak Symbolic and SymbolicWork usage */ |
|
|
688 sym, /* symmetry of pattern */ |
|
|
689 dnum_mem_init_usage, /* min Numeric->Memory for UMF_kernel_init */ |
|
|
690 amd_lunz, /* nz in LU for AMD, with symmetric pivoting */ |
|
|
691 lunz_bound ; /* max nx in LU, for arbitrary row pivoting */ |
|
|
692 |
|
|
693 Int valid, /* set to SYMBOLIC_VALID, for validity check */ |
|
|
694 max_nchains, |
|
|
695 nchains, |
|
|
696 *Chain_start, |
|
|
697 *Chain_maxrows, |
|
|
698 *Chain_maxcols, |
|
|
699 maxnrows, /* largest number of rows in any front */ |
|
|
700 maxncols, /* largest number of columns in any front */ |
|
|
701 *Front_npivcol, /* Front_npivcol [j] = size of jth supercolumn*/ |
|
|
702 *Front_1strow, /* first row index in front j */ |
|
|
703 *Front_leftmostdesc, /* leftmost desc of front j */ |
|
|
704 *Front_parent, /* super-column elimination tree */ |
|
|
705 *Cperm_init, /* initial column ordering */ |
|
|
706 *Rperm_init, /* initial row ordering */ |
|
|
707 *Cdeg, *Rdeg, |
|
|
708 *Esize, |
|
|
709 dense_row_threshold, |
|
|
710 n1, /* number of singletons */ |
|
|
711 nempty, /* MIN (nempty_row, nempty_col) */ |
|
|
712 *Diagonal_map, /* initial "diagonal" (after 2by2) */ |
|
|
713 esize, /* size of Esize array */ |
|
|
714 nfr, |
|
|
715 n_row, n_col, /* matrix A is n_row-by-n_col */ |
|
|
716 nz, /* nz of original matrix */ |
|
|
717 nb, /* block size for BLAS 3 */ |
|
|
718 num_mem_init_usage, /* min Numeric->Memory for UMF_kernel_init */ |
|
|
719 nempty_row, nempty_col, |
|
|
720 |
|
|
721 strategy, |
|
|
722 ordering, |
|
|
723 fixQ, |
|
|
724 prefer_diagonal, |
|
|
725 nzaat, |
|
|
726 nzdiag, |
|
|
727 amd_dmax ; |
|
|
728 |
|
|
729 } SymbolicType ; |
|
|
730 |
|
|
731 |
|
|
732 /* -------------------------------------------------------------------------- */ |
|
|
733 /* for debugging only: */ |
|
|
734 /* -------------------------------------------------------------------------- */ |
|
|
735 |
|
|
736 #include "umf_dump.h" |
|
|
737 |
|
|
738 /* -------------------------------------------------------------------------- */ |
|
|
739 /* for statement coverage testing only: */ |
|
|
740 /* -------------------------------------------------------------------------- */ |
|
|
741 |
|
|
742 #ifdef TESTING |
|
|
743 |
|
|
744 /* for testing integer overflow: */ |
|
|
745 #ifdef TEST_FOR_INTEGER_OVERFLOW |
|
|
746 #undef MAX_MARK |
|
|
747 #define MAX_MARK(n) (3*(n)) |
|
|
748 #endif |
|
|
749 |
|
|
750 /* for testing out-of-memory conditions: */ |
|
|
751 #define UMF_TCOV_TEST |
|
|
752 GLOBAL extern int umf_fail, umf_fail_lo, umf_fail_hi ; |
|
|
753 GLOBAL extern int umf_realloc_fail, umf_realloc_lo, umf_realloc_hi ; |
|
|
754 |
|
|
755 /* for testing malloc count: */ |
|
|
756 #define UMF_MALLOC_COUNT |
|
|
757 |
|
|
758 #endif |
|
|
759 |
|
|
760 #endif |
