|
5164
|
1 /* ========================================================================= */ |
|
|
2 /* === AMD: approximate minimum degree ordering =========================== */ |
|
|
3 /* ========================================================================= */ |
|
|
4 |
|
|
5 /* ------------------------------------------------------------------------- */ |
|
|
6 /* AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis, */ |
|
|
7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README for License. */ |
|
|
8 /* email: davis@cise.ufl.edu CISE Department, Univ. of Florida. */ |
|
|
9 /* web: http://www.cise.ufl.edu/research/sparse/amd */ |
|
|
10 /* ------------------------------------------------------------------------- */ |
|
|
11 |
|
|
12 /* AMD finds a symmetric ordering P of a matrix A so that the Cholesky |
|
|
13 * factorization of P*A*P' has fewer nonzeros and takes less work than the |
|
|
14 * Cholesky factorization of A. If A is not symmetric, then it performs its |
|
|
15 * ordering on the matrix A+A'. Two sets of user-callable routines are |
|
|
16 * provided, one for "int" integers and the other for "long" integers. |
|
|
17 * |
|
|
18 * The method is based on the approximate minimum degree algorithm, discussed |
|
|
19 * in Amestoy, Davis, and Duff, "An approximate degree ordering algorithm", |
|
|
20 * SIAM Journal of Matrix Analysis and Applications, vol. 17, no. 4, pp. |
|
|
21 * 886-905, 1996. This package can perform both the AMD ordering (with |
|
|
22 * aggressive absorption), and the AMDBAR ordering (without aggressive |
|
|
23 * absorption) discussed in the above paper. This package differs from the |
|
|
24 * Fortran codes discussed in the paper: |
|
|
25 * |
|
|
26 * (1) it can ignore "dense" rows and columns, leading to faster run times |
|
|
27 * (2) it computes the ordering of A+A' if A is not symmetric |
|
|
28 * (3) it is followed by a depth-first post-ordering of the assembly tree |
|
|
29 * (or supernodal elimination tree) |
|
|
30 * |
|
|
31 * For historical reasons, the Fortran versions, amd.f and amdbar.f, have |
|
|
32 * been left (nearly) unchanged. They compute the identical ordering as |
|
|
33 * described in the above paper. |
|
|
34 */ |
|
|
35 |
|
|
36 #ifndef AMD_H |
|
|
37 #define AMD_H |
|
|
38 |
|
|
39 int amd_order ( /* returns 0 if OK, negative value if error */ |
|
|
40 int n, /* A is n-by-n. n must be >= 0. */ |
|
|
41 const int Ap [ ], /* column pointers for A, of size n+1 */ |
|
|
42 const int Ai [ ], /* row indices of A, of size nz = Ap [n] */ |
|
|
43 int P [ ], /* output permutation, of size n */ |
|
|
44 double Control [ ], /* input Control settings, of size AMD_CONTROL */ |
|
|
45 double Info [ ] /* output Info statistics, of size AMD_INFO */ |
|
|
46 ) ; |
|
|
47 |
|
|
48 long amd_l_order ( /* see above for description of arguments */ |
|
|
49 long n, |
|
|
50 const long Ap [ ], |
|
|
51 const long Ai [ ], |
|
|
52 long P [ ], |
|
|
53 double Control [ ], |
|
|
54 double Info [ ] |
|
|
55 ) ; |
|
|
56 |
|
|
57 /* Input arguments (not modified): |
|
|
58 * |
|
|
59 * n: the matrix A is n-by-n. |
|
|
60 * Ap: an int/long array of size n+1, containing the column pointers of A. |
|
|
61 * Ai: an int/long array of size nz, containing the row indices of A, |
|
|
62 * where nz = Ap [n]. |
|
|
63 * Control: a double array of size AMD_CONTROL, containing control |
|
|
64 * parameters. Defaults are used if Control is NULL. |
|
|
65 * |
|
|
66 * Output arguments (not defined on input): |
|
|
67 * |
|
|
68 * P: an int/long array of size n, containing the output permutation. If |
|
|
69 * row i is the kth pivot row, then P [k] = i. In MATLAB notation, |
|
|
70 * the reordered matrix is A (P,P). |
|
|
71 * Info: a double array of size AMD_INFO, containing statistical |
|
|
72 * information. Ignored if Info is NULL. |
|
|
73 * |
|
|
74 * On input, the matrix A is stored in column-oriented form. The row indices |
|
|
75 * of nonzero entries in column j are stored in Ai [Ap [j] ... Ap [j+1]-1]. |
|
|
76 * The row indices must appear in ascending order in each column, and there |
|
|
77 * must not be any duplicate entries. Row indices must be in the range 0 to |
|
|
78 * n-1. Ap [0] must be zero, and thus nz = Ap [n] is the number of nonzeros |
|
|
79 * in A. The array Ap is of size n+1, and the array Ai is of size nz = Ap [n]. |
|
|
80 * The matrix does not need to be symmetric, and the diagonal does not need to |
|
|
81 * be present (if diagonal entries are present, they are ignored except for |
|
|
82 * the output statistic Info [AMD_NZDIAG]). The arrays Ai and Ap are not |
|
|
83 * modified. This form of the Ap and Ai arrays to represent the nonzero |
|
|
84 * pattern of the matrix A is the same as that used internally by MATLAB. |
|
|
85 * If you wish to use a more flexible input structure, please see the |
|
|
86 * umfpack_*_triplet_to_col routines in the UMFPACK package, at |
|
|
87 * http://www.cise.ufl.edu/research/sparse/umfpack, or use the amd_preprocess |
|
|
88 * routine discussed below. |
|
|
89 * |
|
|
90 * Restrictions: n >= 0. Ap [0] = 0. Ap [j] <= Ap [j+1] for all j in the |
|
|
91 * range 0 to n-1. nz = Ap [n] >= 0. For all j in the range 0 to n-1, |
|
|
92 * and for all p in the range Ap [j] to Ap [j+1]-2, Ai [p] < Ai [p+1] must |
|
|
93 * hold. Ai [0..nz-1] must be in the range 0 to n-1. To avoid integer |
|
|
94 * overflow, (2.4*nz + 8*n) < INT_MAX / sizeof (int) for must hold for the |
|
|
95 * "int" version. (2.4*nz + 8*n) < LONG_MAX / sizeof (long) must hold |
|
|
96 * for the "long" version. Finally, Ai, Ap, and P must not be NULL. If |
|
|
97 * any of these restrictions are not met, AMD returns AMD_INVALID. |
|
|
98 * |
|
|
99 * AMD returns: |
|
|
100 * |
|
|
101 * AMD_OK if the matrix is valid and sufficient memory can be allocated to |
|
|
102 * perform the ordering. |
|
|
103 * |
|
|
104 * AMD_OUT_OF_MEMORY if not enough memory can be allocated. |
|
|
105 * |
|
|
106 * AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if P is |
|
|
107 * NULL. |
|
|
108 * |
|
|
109 * The AMD routine first forms the pattern of the matrix A+A', and then |
|
|
110 * computes a fill-reducing ordering, P. If P [k] = i, then row/column i of |
|
|
111 * the original is the kth pivotal row. In MATLAB notation, the permuted |
|
|
112 * matrix is A (P,P), except that 0-based indexing is used instead of the |
|
|
113 * 1-based indexing in MATLAB. |
|
|
114 * |
|
|
115 * The Control array is used to set various parameters for AMD. If a NULL |
|
|
116 * pointer is passed, default values are used. The Control array is not |
|
|
117 * modified. |
|
|
118 * |
|
|
119 * Control [AMD_DENSE]: controls the threshold for "dense" rows/columns. |
|
|
120 * A dense row/column in A+A' can cause AMD to spend a lot of time in |
|
|
121 * ordering the matrix. If Control [AMD_DENSE] >= 0, rows/columns |
|
|
122 * with more than Control [AMD_DENSE] * sqrt (n) entries are ignored |
|
|
123 * during the ordering, and placed last in the output order. The |
|
|
124 * default value of Control [AMD_DENSE] is 10. If negative, no |
|
|
125 * rows/columns are treated as "dense". Rows/columns with 16 or |
|
|
126 * fewer off-diagonal entries are never considered "dense". |
|
|
127 * |
|
|
128 * Control [AMD_AGGRESSIVE]: controls whether or not to use aggressive |
|
|
129 * absorption, in which a prior element is absorbed into the current |
|
|
130 * element if is a subset of the current element, even if it is not |
|
|
131 * adjacent to the current pivot element (refer to Amestoy, Davis, |
|
|
132 * & Duff, 1996, for more details). The default value is nonzero, |
|
|
133 * which means to perform aggressive absorption. This nearly always |
|
|
134 * leads to a better ordering (because the approximate degrees are |
|
|
135 * more accurate) and a lower execution time. There are cases where |
|
|
136 * it can lead to a slightly worse ordering, however. To turn it off, |
|
|
137 * set Control [AMD_AGGRESSIVE] to 0. |
|
|
138 * |
|
|
139 * Control [2..4] are not used in the current version, but may be used in |
|
|
140 * future versions. |
|
|
141 * |
|
|
142 * The Info array provides statistics about the ordering on output. If it is |
|
|
143 * not present, the statistics are not returned. This is not an error |
|
|
144 * condition. |
|
|
145 * |
|
|
146 * Info [AMD_STATUS]: the return value of AMD, either AMD_OK, |
|
|
147 * AMD_OUT_OF_MEMORY, or AMD_INVALID. |
|
|
148 * |
|
|
149 * Info [AMD_N]: n, the size of the input matrix |
|
|
150 * |
|
|
151 * Info [AMD_NZ]: the number of nonzeros in A, nz = Ap [n] |
|
|
152 * |
|
|
153 * Info [AMD_SYMMETRY]: the symmetry of the matrix A. It is the number |
|
|
154 * of "matched" off-diagonal entries divided by the total number of |
|
|
155 * off-diagonal entries. An entry A(i,j) is matched if A(j,i) is also |
|
|
156 * an entry, for any pair (i,j) for which i != j. In MATLAB notation, |
|
|
157 * S = spones (A) ; |
|
|
158 * B = tril (S, -1) + triu (S, 1) ; |
|
|
159 * symmetry = nnz (B & B') / nnz (B) ; |
|
|
160 * |
|
|
161 * Info [AMD_NZDIAG]: the number of entries on the diagonal of A. |
|
|
162 * |
|
|
163 * Info [AMD_NZ_A_PLUS_AT]: the number of nonzeros in A+A', excluding the |
|
|
164 * diagonal. If A is perfectly symmetric (Info [AMD_SYMMETRY] = 1) |
|
|
165 * with a fully nonzero diagonal, then Info [AMD_NZ_A_PLUS_AT] = nz-n |
|
|
166 * (the smallest possible value). If A is perfectly unsymmetric |
|
|
167 * (Info [AMD_SYMMETRY] = 0, for an upper triangular matrix, for |
|
|
168 * example) with no diagonal, then Info [AMD_NZ_A_PLUS_AT] = 2*nz |
|
|
169 * (the largest possible value). |
|
|
170 * |
|
|
171 * Info [AMD_NDENSE]: the number of "dense" rows/columns of A+A' that were |
|
|
172 * removed from A prior to ordering. These are placed last in the |
|
|
173 * output order P. |
|
|
174 * |
|
|
175 * Info [AMD_MEMORY]: the amount of memory used by AMD, in bytes. In the |
|
|
176 * current version, this is 1.2 * Info [AMD_NZ_A_PLUS_AT] + 9*n |
|
|
177 * times the size of an integer. This is at most 2.4nz + 9n. This |
|
|
178 * excludes the size of the input arguments Ai, Ap, and P, which have |
|
|
179 * a total size of nz + 2*n + 1 integers. |
|
|
180 * |
|
|
181 * Info [AMD_NCMPA]: the number of garbage collections performed. |
|
|
182 * |
|
|
183 * Info [AMD_LNZ]: the number of nonzeros in L (excluding the diagonal). |
|
|
184 * This is a slight upper bound because mass elimination is combined |
|
|
185 * with the approximate degree update. It is a rough upper bound if |
|
|
186 * there are many "dense" rows/columns. The rest of the statistics, |
|
|
187 * below, are also slight or rough upper bounds, for the same reasons. |
|
|
188 * The post-ordering of the assembly tree might also not exactly |
|
|
189 * correspond to a true elimination tree postordering. |
|
|
190 * |
|
|
191 * Info [AMD_NDIV]: the number of divide operations for a subsequent LDL' |
|
|
192 * or LU factorization of the permuted matrix A (P,P). |
|
|
193 * |
|
|
194 * Info [AMD_NMULTSUBS_LDL]: the number of multiply-subtract pairs for a |
|
|
195 * subsequent LDL' factorization of A (P,P). |
|
|
196 * |
|
|
197 * Info [AMD_NMULTSUBS_LU]: the number of multiply-subtract pairs for a |
|
|
198 * subsequent LU factorization of A (P,P), assuming that no numerical |
|
|
199 * pivoting is required. |
|
|
200 * |
|
|
201 * Info [AMD_DMAX]: the maximum number of nonzeros in any column of L, |
|
|
202 * including the diagonal. |
|
|
203 * |
|
|
204 * Info [14..19] are not used in the current version, but may be used in |
|
|
205 * future versions. |
|
|
206 */ |
|
|
207 |
|
|
208 /* ------------------------------------------------------------------------- */ |
|
|
209 /* AMD preprocess */ |
|
|
210 /* ------------------------------------------------------------------------- */ |
|
|
211 |
|
|
212 /* amd_preprocess: sorts, removes duplicate entries, and transposes the |
|
|
213 * nonzero pattern of a column-form matrix A, to obtain the matrix R. |
|
|
214 * |
|
|
215 * Alternatively, you can consider this routine as constructing a row-form |
|
|
216 * matrix from a column-form matrix. Duplicate entries are allowed in A (and |
|
|
217 * removed in R). The columns of R are sorted. Checks its input A for errors. |
|
|
218 * |
|
|
219 * On input, A can have unsorted columns, and can have duplicate entries. |
|
|
220 * Ap [0] must still be zero, and Ap must be monotonically nondecreasing. |
|
|
221 * Row indices must be in the range 0 to n-1. |
|
|
222 * |
|
|
223 * On output, if this routine returns AMD_OK, then the matrix R is a valid |
|
|
224 * input matrix for AMD_order. It has sorted columns, with no duplicate |
|
|
225 * entries in each column. Since AMD_order operates on the matrix A+A', it |
|
|
226 * can just as easily use A or A', so the transpose has no significant effect |
|
|
227 * (except for minor tie-breaking, which can lead to a minor effect in the |
|
|
228 * quality of the ordering). As an example, compare the output of amd_demo.c |
|
|
229 * and amd_demo2.c. |
|
|
230 * |
|
|
231 * This routine transposes A to get R because that's the simplest way to |
|
|
232 * sort and remove duplicate entries from a matrix. |
|
|
233 * |
|
|
234 * Allocates 2*n integer work arrays, and free's them when done. |
|
|
235 * |
|
|
236 * If you wish to call amd_order, but do not know if your matrix has unsorted |
|
|
237 * columns or duplicate entries, then you can use the following code, which is |
|
|
238 * fairly efficient. amd_order will not allocate any internal matrix until |
|
|
239 * it checks that the input matrix is valid, so the method below is memory- |
|
|
240 * efficient as well. This code snippet assumes that Rp and Ri are already |
|
|
241 * allocated, and are the same size as Ap and Ai respectively. |
|
|
242 |
|
|
243 result = amd_order (n, p, Ap, Ai, Control, Info) ; |
|
|
244 if (result == AMD_INVALID) |
|
|
245 { |
|
|
246 if (amd_preprocess (n, Ap, Ai, Rp, Ri) == AMD_OK) |
|
|
247 { |
|
|
248 result = amd_order (n, p, Rp, Ri, Control, Info) ; |
|
|
249 } |
|
|
250 } |
|
|
251 |
|
|
252 * amd_preprocess will still return AMD_INVALID if any row index in Ai is out |
|
|
253 * of range or if the Ap array is invalid. These errors are not corrected by |
|
|
254 * amd_preprocess since they represent a more serious error that should be |
|
|
255 * flagged with the AMD_INVALID error code. |
|
|
256 */ |
|
|
257 |
|
|
258 int amd_preprocess |
|
|
259 ( |
|
|
260 int n, |
|
|
261 const int Ap [ ], |
|
|
262 const int Ai [ ], |
|
|
263 int Rp [ ], |
|
|
264 int Ri [ ] |
|
|
265 ) ; |
|
|
266 |
|
|
267 long amd_l_preprocess |
|
|
268 ( |
|
|
269 long n, |
|
|
270 const long Ap [ ], |
|
|
271 const long Ai [ ], |
|
|
272 long Rp [ ], |
|
|
273 long Ri [ ] |
|
|
274 ) ; |
|
|
275 |
|
|
276 /* Input arguments (not modified): |
|
|
277 * |
|
|
278 * n: the matrix A is n-by-n. |
|
|
279 * Ap: an int/long array of size n+1, containing the column pointers of A. |
|
|
280 * Ai: an int/long array of size nz, containing the row indices of A, |
|
|
281 * where nz = Ap [n]. |
|
|
282 * The nonzero pattern of column j of A is in Ai [Ap [j] ... Ap [j+1]-1]. |
|
|
283 * Ap [0] must be zero, and Ap [j] <= Ap [j+1] must hold for all j in the |
|
|
284 * range 0 to n-1. Row indices in Ai must be in the range 0 to n-1. |
|
|
285 * The row indices in any one column need not be sorted, and duplicates |
|
|
286 * may exist. |
|
|
287 * |
|
|
288 * Output arguments (not defined on input): |
|
|
289 * |
|
|
290 * Rp: an int/long array of size n+1, containing the column pointers of R. |
|
|
291 * Ri: an int/long array of size rnz, containing the row indices of R, |
|
|
292 * where rnz = Rp [n]. Note that Rp [n] will be less than Ap [n] if |
|
|
293 * duplicates appear in A. In general, Rp [n] <= Ap [n]. |
|
|
294 * The data structure for R is the same as A, except that each column of |
|
|
295 * R contains sorted row indices, and no duplicates appear in any column. |
|
|
296 * |
|
|
297 * amd_preprocess returns: |
|
|
298 * |
|
|
299 * AMD_OK if the matrix A is valid and sufficient memory can be allocated |
|
|
300 * to perform the preprocessing. |
|
|
301 * |
|
|
302 * AMD_OUT_OF_MEMORY if not enough memory can be allocated. |
|
|
303 * |
|
|
304 * AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if Rp or |
|
|
305 * Ri are NULL. |
|
|
306 */ |
|
|
307 |
|
|
308 /* ------------------------------------------------------------------------- */ |
|
|
309 /* AMD Control and Info arrays */ |
|
|
310 /* ------------------------------------------------------------------------- */ |
|
|
311 |
|
|
312 /* amd_defaults: sets the default control settings */ |
|
|
313 void amd_defaults (double Control [ ]) ; |
|
|
314 void amd_l_defaults (double Control [ ]) ; |
|
|
315 |
|
|
316 /* amd_control: prints the control settings */ |
|
|
317 void amd_control (double Control [ ]) ; |
|
|
318 void amd_l_control (double Control [ ]) ; |
|
|
319 |
|
|
320 /* amd_info: prints the statistics */ |
|
|
321 void amd_info (double Info [ ]) ; |
|
|
322 void amd_l_info (double Info [ ]) ; |
|
|
323 |
|
|
324 #define AMD_CONTROL 5 /* size of Control array */ |
|
|
325 #define AMD_INFO 20 /* size of Info array */ |
|
|
326 |
|
|
327 /* contents of Control */ |
|
|
328 #define AMD_DENSE 0 /* "dense" if degree > Control [0] * sqrt (n) */ |
|
|
329 #define AMD_AGGRESSIVE 1 /* do aggressive absorption if Control [1] != 0 */ |
|
|
330 |
|
|
331 /* default Control settings */ |
|
|
332 #define AMD_DEFAULT_DENSE 10.0 /* default "dense" degree 10*sqrt(n) */ |
|
|
333 #define AMD_DEFAULT_AGGRESSIVE 1 /* do aggressive absorption by default */ |
|
|
334 |
|
|
335 /* contents of Info */ |
|
|
336 #define AMD_STATUS 0 /* return value of amd_order and amd_l_order */ |
|
|
337 #define AMD_N 1 /* A is n-by-n */ |
|
|
338 #define AMD_NZ 2 /* number of nonzeros in A */ |
|
|
339 #define AMD_SYMMETRY 3 /* symmetry of pattern (1 is sym., 0 is unsym.) */ |
|
|
340 #define AMD_NZDIAG 4 /* # of entries on diagonal */ |
|
|
341 #define AMD_NZ_A_PLUS_AT 5 /* nz in A+A' */ |
|
|
342 #define AMD_NDENSE 6 /* number of "dense" rows/columns in A */ |
|
|
343 #define AMD_MEMORY 7 /* amount of memory used by AMD */ |
|
|
344 #define AMD_NCMPA 8 /* number of garbage collections in AMD */ |
|
|
345 #define AMD_LNZ 9 /* approx. nz in L, excluding the diagonal */ |
|
|
346 #define AMD_NDIV 10 /* number of fl. point divides for LU and LDL' */ |
|
|
347 #define AMD_NMULTSUBS_LDL 11 /* number of fl. point (*,-) pairs for LDL' */ |
|
|
348 #define AMD_NMULTSUBS_LU 12 /* number of fl. point (*,-) pairs for LU */ |
|
|
349 #define AMD_DMAX 13 /* max nz. in any column of L, incl. diagonal */ |
|
|
350 |
|
|
351 /* ------------------------------------------------------------------------- */ |
|
|
352 /* return values of AMD */ |
|
|
353 /* ------------------------------------------------------------------------- */ |
|
|
354 |
|
|
355 #define AMD_OK 0 /* success */ |
|
|
356 #define AMD_OUT_OF_MEMORY -1 /* malloc failed */ |
|
|
357 #define AMD_INVALID -2 /* input arguments are not valid */ |
|
|
358 |
|
|
359 #endif |
