1 /* ========================================================================= */ 2 /* === AMD: approximate minimum degree ordering =========================== */ 3 /* ========================================================================= */ 4 5 /* ------------------------------------------------------------------------- */ 6 /* AMD Version 2.2, Copyright (c) 2007 by Timothy A. Davis, */ 7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */ 8 /* email: DrTimothyAldenDavis@gmail.com */ 9 /* ------------------------------------------------------------------------- */ 10 11 /* AMD finds a symmetric ordering P of a matrix A so that the Cholesky 12 * factorization of P*A*P' has fewer nonzeros and takes less work than the 13 * Cholesky factorization of A. If A is not symmetric, then it performs its 14 * ordering on the matrix A+A'. Two sets of user-callable routines are 15 * provided, one for int integers and the other for SuiteSparse_long integers. 16 * 17 * The method is based on the approximate minimum degree algorithm, discussed 18 * in Amestoy, Davis, and Duff, "An approximate degree ordering algorithm", 19 * SIAM Journal of Matrix Analysis and Applications, vol. 17, no. 4, pp. 20 * 886-905, 1996. This package can perform both the AMD ordering (with 21 * aggressive absorption), and the AMDBAR ordering (without aggressive 22 * absorption) discussed in the above paper. This package differs from the 23 * Fortran codes discussed in the paper: 24 * 25 * (1) it can ignore "dense" rows and columns, leading to faster run times 26 * (2) it computes the ordering of A+A' if A is not symmetric test03()27 * (3) it is followed by a depth-first post-ordering of the assembly tree 28 * (or supernodal elimination tree) 29 * 30 * For historical reasons, the Fortran versions, amd.f and amdbar.f, have 31 * been left (nearly) unchanged. They compute the identical ordering as 32 * described in the above paper. 33 */ 34 35 #ifndef AMD_H 36 #define AMD_H 37 38 /* make it easy for C++ programs to include AMD */ 39 #ifdef __cplusplus 40 extern "C" { 41 #endif 42 43 /* get the definition of size_t: */ 44 #include <stddef.h> 45 46 #include "SuiteSparse_config.h" 47 48 int amd_order /* returns AMD_OK, AMD_OK_BUT_JUMBLED, 49 * AMD_INVALID, or AMD_OUT_OF_MEMORY */ 50 ( 51 int n, /* A is n-by-n. n must be >= 0. */ 52 const int Ap [ ], /* column pointers for A, of size n+1 */ 53 const int Ai [ ], /* row indices of A, of size nz = Ap [n] */ 54 int P [ ], /* output permutation, of size n */ 55 double Control [ ], /* input Control settings, of size AMD_CONTROL */ 56 double Info [ ] /* output Info statistics, of size AMD_INFO */ 57 ) ; 58 59 SuiteSparse_long amd_l_order /* see above for description of arguments */ 60 ( 61 SuiteSparse_long n, 62 const SuiteSparse_long Ap [ ], 63 const SuiteSparse_long Ai [ ], 64 SuiteSparse_long P [ ], 65 double Control [ ], 66 double Info [ ] 67 ) ; 68 69 /* Input arguments (not modified): 70 * 71 * n: the matrix A is n-by-n. 72 * Ap: an int/SuiteSparse_long array of size n+1, containing column 73 * pointers of A. 74 * Ai: an int/SuiteSparse_long array of size nz, containing the row 75 * indices of A, where nz = Ap [n]. 76 * Control: a double array of size AMD_CONTROL, containing control 77 * parameters. Defaults are used if Control is NULL. 78 * 79 * Output arguments (not defined on input): 80 * 81 * P: an int/SuiteSparse_long array of size n, containing the output 82 * permutation. If row i is the kth pivot row, then P [k] = i. In 83 * MATLAB notation, the reordered matrix is A (P,P). 84 * Info: a double array of size AMD_INFO, containing statistical 85 * information. Ignored if Info is NULL. 86 * 87 * On input, the matrix A is stored in column-oriented form. The row indices 88 * of nonzero entries in column j are stored in Ai [Ap [j] ... Ap [j+1]-1]. 89 * 90 * If the row indices appear in ascending order in each column, and there 91 * are no duplicate entries, then amd_order is slightly more efficient in 92 * terms of time and memory usage. If this condition does not hold, a copy 93 * of the matrix is created (where these conditions do hold), and the copy is 94 * ordered. This feature is new to v2.0 (v1.2 and earlier required this 95 * condition to hold for the input matrix). 96 * 97 * Row indices must be in the range 0 to 98 * n-1. Ap [0] must be zero, and thus nz = Ap [n] is the number of nonzeros 99 * in A. The array Ap is of size n+1, and the array Ai is of size nz = Ap [n]. 100 * The matrix does not need to be symmetric, and the diagonal does not need to 101 * be present (if diagonal entries are present, they are ignored except for 102 * the output statistic Info [AMD_NZDIAG]). The arrays Ai and Ap are not 103 * modified. This form of the Ap and Ai arrays to represent the nonzero 104 * pattern of the matrix A is the same as that used internally by MATLAB. 105 * If you wish to use a more flexible input structure, please see the 106 * umfpack_*_triplet_to_col routines in the UMFPACK package, at 107 * http://www.suitesparse.com. 108 * 109 * Restrictions: n >= 0. Ap [0] = 0. Ap [j] <= Ap [j+1] for all j in the 110 * range 0 to n-1. nz = Ap [n] >= 0. Ai [0..nz-1] must be in the range 0 111 * to n-1. Finally, Ai, Ap, and P must not be NULL. If any of these 112 * restrictions are not met, AMD returns AMD_INVALID. 113 * 114 * AMD returns: 115 * 116 * AMD_OK if the matrix is valid and sufficient memory can be allocated to 117 * perform the ordering. 118 * 119 * AMD_OUT_OF_MEMORY if not enough memory can be allocated. 120 * 121 * AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if P is 122 * NULL. 123 * 124 * AMD_OK_BUT_JUMBLED if the matrix had unsorted columns, and/or duplicate 125 * entries, but was otherwise valid. 126 * 127 * The AMD routine first forms the pattern of the matrix A+A', and then 128 * computes a fill-reducing ordering, P. If P [k] = i, then row/column i of 129 * the original is the kth pivotal row. In MATLAB notation, the permuted 130 * matrix is A (P,P), except that 0-based indexing is used instead of the 131 * 1-based indexing in MATLAB. 132 * 133 * The Control array is used to set various parameters for AMD. If a NULL 134 * pointer is passed, default values are used. The Control array is not 135 * modified. 136 * 137 * Control [AMD_DENSE]: controls the threshold for "dense" rows/columns. 138 * A dense row/column in A+A' can cause AMD to spend a lot of time in 139 * ordering the matrix. If Control [AMD_DENSE] >= 0, rows/columns 140 * with more than Control [AMD_DENSE] * sqrt (n) entries are ignored 141 * during the ordering, and placed last in the output order. The 142 * default value of Control [AMD_DENSE] is 10. If negative, no 143 * rows/columns are treated as "dense". Rows/columns with 16 or 144 * fewer off-diagonal entries are never considered "dense". 145 * 146 * Control [AMD_AGGRESSIVE]: controls whether or not to use aggressive 147 * absorption, in which a prior element is absorbed into the current 148 * element if is a subset of the current element, even if it is not 149 * adjacent to the current pivot element (refer to Amestoy, Davis, 150 * & Duff, 1996, for more details). The default value is nonzero, 151 * which means to perform aggressive absorption. This nearly always 152 * leads to a better ordering (because the approximate degrees are 153 * more accurate) and a lower execution time. There are cases where 154 * it can lead to a slightly worse ordering, however. To turn it off, 155 * set Control [AMD_AGGRESSIVE] to 0. 156 * 157 * Control [2..4] are not used in the current version, but may be used in 158 * future versions. 159 * 160 * The Info array provides statistics about the ordering on output. If it is 161 * not present, the statistics are not returned. This is not an error 162 * condition. 163 * 164 * Info [AMD_STATUS]: the return value of AMD, either AMD_OK, 165 * AMD_OK_BUT_JUMBLED, AMD_OUT_OF_MEMORY, or AMD_INVALID. 166 * 167 * Info [AMD_N]: n, the size of the input matrix 168 * 169 * Info [AMD_NZ]: the number of nonzeros in A, nz = Ap [n] 170 * 171 * Info [AMD_SYMMETRY]: the symmetry of the matrix A. It is the number 172 * of "matched" off-diagonal entries divided by the total number of 173 * off-diagonal entries. An entry A(i,j) is matched if A(j,i) is also 174 * an entry, for any pair (i,j) for which i != j. In MATLAB notation, 175 * S = spones (A) ; 176 * B = tril (S, -1) + triu (S, 1) ; 177 * symmetry = nnz (B & B') / nnz (B) ; 178 * 179 * Info [AMD_NZDIAG]: the number of entries on the diagonal of A. 180 * 181 * Info [AMD_NZ_A_PLUS_AT]: the number of nonzeros in A+A', excluding the 182 * diagonal. If A is perfectly symmetric (Info [AMD_SYMMETRY] = 1) 183 * with a fully nonzero diagonal, then Info [AMD_NZ_A_PLUS_AT] = nz-n 184 * (the smallest possible value). If A is perfectly unsymmetric 185 * (Info [AMD_SYMMETRY] = 0, for an upper triangular matrix, for 186 * example) with no diagonal, then Info [AMD_NZ_A_PLUS_AT] = 2*nz 187 * (the largest possible value). 188 * 189 * Info [AMD_NDENSE]: the number of "dense" rows/columns of A+A' that were 190 * removed from A prior to ordering. These are placed last in the 191 * output order P. 192 * 193 * Info [AMD_MEMORY]: the amount of memory used by AMD, in bytes. In the 194 * current version, this is 1.2 * Info [AMD_NZ_A_PLUS_AT] + 9*n 195 * times the size of an integer. This is at most 2.4nz + 9n. This 196 * excludes the size of the input arguments Ai, Ap, and P, which have 197 * a total size of nz + 2*n + 1 integers. 198 * 199 * Info [AMD_NCMPA]: the number of garbage collections performed. 200 * 201 * Info [AMD_LNZ]: the number of nonzeros in L (excluding the diagonal). 202 * This is a slight upper bound because mass elimination is combined 203 * with the approximate degree update. It is a rough upper bound if 204 * there are many "dense" rows/columns. The rest of the statistics, 205 * below, are also slight or rough upper bounds, for the same reasons. 206 * The post-ordering of the assembly tree might also not exactly 207 * correspond to a true elimination tree postordering. 208 * 209 * Info [AMD_NDIV]: the number of divide operations for a subsequent LDL' 210 * or LU factorization of the permuted matrix A (P,P). 211 * 212 * Info [AMD_NMULTSUBS_LDL]: the number of multiply-subtract pairs for a 213 * subsequent LDL' factorization of A (P,P). 214 * 215 * Info [AMD_NMULTSUBS_LU]: the number of multiply-subtract pairs for a 216 * subsequent LU factorization of A (P,P), assuming that no numerical 217 * pivoting is required. 218 * 219 * Info [AMD_DMAX]: the maximum number of nonzeros in any column of L, 220 * including the diagonal. 221 * 222 * Info [14..19] are not used in the current version, but may be used in 223 * future versions. 224 */ 225 226 /* ------------------------------------------------------------------------- */ 227 /* direct interface to AMD */ 228 /* ------------------------------------------------------------------------- */ 229 230 /* amd_2 is the primary AMD ordering routine. It is not meant to be 231 * user-callable because of its restrictive inputs and because it destroys 232 * the user's input matrix. It does not check its inputs for errors, either. 233 * However, if you can work with these restrictions it can be faster than 234 * amd_order and use less memory (assuming that you can create your own copy 235 * of the matrix for AMD to destroy). Refer to AMD/Source/amd_2.c for a 236 * description of each parameter. */ 237 238 void amd_2 239 ( 240 int n, 241 int Pe [ ], 242 int Iw [ ], 243 int Len [ ], 244 int iwlen, 245 int pfree, 246 int Nv [ ], 247 int Next [ ], 248 int Last [ ], 249 int Head [ ], 250 int Elen [ ], 251 int Degree [ ], 252 int W [ ], 253 double Control [ ], 254 double Info [ ] 255 ) ; 256 257 void amd_l2 258 ( 259 SuiteSparse_long n, 260 SuiteSparse_long Pe [ ], 261 SuiteSparse_long Iw [ ], 262 SuiteSparse_long Len [ ], 263 SuiteSparse_long iwlen, 264 SuiteSparse_long pfree, 265 SuiteSparse_long Nv [ ], 266 SuiteSparse_long Next [ ], 267 SuiteSparse_long Last [ ], 268 SuiteSparse_long Head [ ], 269 SuiteSparse_long Elen [ ], 270 SuiteSparse_long Degree [ ], 271 SuiteSparse_long W [ ], 272 double Control [ ], 273 double Info [ ] 274 ) ; 275 276 /* ------------------------------------------------------------------------- */ 277 /* amd_valid */ 278 /* ------------------------------------------------------------------------- */ 279 280 /* Returns AMD_OK or AMD_OK_BUT_JUMBLED if the matrix is valid as input to 281 * amd_order; the latter is returned if the matrix has unsorted and/or 282 * duplicate row indices in one or more columns. Returns AMD_INVALID if the 283 * matrix cannot be passed to amd_order. For amd_order, the matrix must also 284 * be square. The first two arguments are the number of rows and the number 285 * of columns of the matrix. For its use in AMD, these must both equal n. 286 * 287 * NOTE: this routine returned TRUE/FALSE in v1.2 and earlier. 288 */ 289 290 int amd_valid 291 ( 292 int n_row, /* # of rows */ 293 int n_col, /* # of columns */ 294 const int Ap [ ], /* column pointers, of size n_col+1 */ 295 const int Ai [ ] /* row indices, of size Ap [n_col] */ 296 ) ; 297 298 SuiteSparse_long amd_l_valid 299 ( 300 SuiteSparse_long n_row, 301 SuiteSparse_long n_col, 302 const SuiteSparse_long Ap [ ], 303 const SuiteSparse_long Ai [ ] 304 ) ; 305 306 /* ------------------------------------------------------------------------- */ 307 /* AMD memory manager and printf routines */ 308 /* ------------------------------------------------------------------------- */ 309 310 /* The user can redefine these to change the malloc, free, and printf routines 311 * that AMD uses. */ 312 313 #ifndef EXTERN 314 #define EXTERN extern 315 #endif 316 317 EXTERN void *(*amd_malloc) (size_t) ; /* pointer to malloc */ 318 EXTERN void (*amd_free) (void *) ; /* pointer to free */ 319 EXTERN void *(*amd_realloc) (void *, size_t) ; /* pointer to realloc */ 320 EXTERN void *(*amd_calloc) (size_t, size_t) ; /* pointer to calloc */ 321 EXTERN int (*amd_printf) (const char *, ...) ; /* pointer to printf */ 322 323 /* ------------------------------------------------------------------------- */ 324 /* AMD Control and Info arrays */ 325 /* ------------------------------------------------------------------------- */ 326 327 /* amd_defaults: sets the default control settings */ 328 void amd_defaults (double Control [ ]) ; 329 void amd_l_defaults (double Control [ ]) ; 330 331 /* amd_control: prints the control settings */ 332 void amd_control (double Control [ ]) ; 333 void amd_l_control (double Control [ ]) ; 334 335 /* amd_info: prints the statistics */ 336 void amd_info (double Info [ ]) ; 337 void amd_l_info (double Info [ ]) ; 338 339 #define AMD_CONTROL 5 /* size of Control array */ 340 #define AMD_INFO 20 /* size of Info array */ 341 342 /* contents of Control */ 343 #define AMD_DENSE 0 /* "dense" if degree > Control [0] * sqrt (n) */ 344 #define AMD_AGGRESSIVE 1 /* do aggressive absorption if Control [1] != 0 */ 345 346 /* default Control settings */ 347 #define AMD_DEFAULT_DENSE 10.0 /* default "dense" degree 10*sqrt(n) */ 348 #define AMD_DEFAULT_AGGRESSIVE 1 /* do aggressive absorption by default */ 349 350 /* contents of Info */ 351 #define AMD_STATUS 0 /* return value of amd_order and amd_l_order */ 352 #define AMD_N 1 /* A is n-by-n */ 353 #define AMD_NZ 2 /* number of nonzeros in A */ 354 #define AMD_SYMMETRY 3 /* symmetry of pattern (1 is sym., 0 is unsym.) */ 355 #define AMD_NZDIAG 4 /* # of entries on diagonal */ 356 #define AMD_NZ_A_PLUS_AT 5 /* nz in A+A' */ 357 #define AMD_NDENSE 6 /* number of "dense" rows/columns in A */ 358 #define AMD_MEMORY 7 /* amount of memory used by AMD */ 359 #define AMD_NCMPA 8 /* number of garbage collections in AMD */ 360 #define AMD_LNZ 9 /* approx. nz in L, excluding the diagonal */ 361 #define AMD_NDIV 10 /* number of fl. point divides for LU and LDL' */ 362 #define AMD_NMULTSUBS_LDL 11 /* number of fl. point (*,-) pairs for LDL' */ 363 #define AMD_NMULTSUBS_LU 12 /* number of fl. point (*,-) pairs for LU */ 364 #define AMD_DMAX 13 /* max nz. in any column of L, incl. diagonal */ 365 366 /* ------------------------------------------------------------------------- */ 367 /* return values of AMD */ 368 /* ------------------------------------------------------------------------- */ 369 370 #define AMD_OK 0 /* success */ 371 #define AMD_OUT_OF_MEMORY -1 /* malloc failed, or problem too large */ 372 #define AMD_INVALID -2 /* input arguments are not valid */ 373 #define AMD_OK_BUT_JUMBLED 1 /* input matrix is OK for amd_order, but 374 * columns were not sorted, and/or duplicate entries were present. AMD had 375 * to do extra work before ordering the matrix. This is a warning, not an 376 * error. */ 377 378 /* ========================================================================== */ 379 /* === AMD version ========================================================== */ 380 /* ========================================================================== */ 381 382 /* AMD Version 1.2 and later include the following definitions. 383 * As an example, to test if the version you are using is 1.2 or later: 384 * 385 * #ifdef AMD_VERSION 386 * if (AMD_VERSION >= AMD_VERSION_CODE (1,2)) ... 387 * #endif 388 * 389 * This also works during compile-time: 390 * 391 * #if defined(AMD_VERSION) && (AMD_VERSION >= AMD_VERSION_CODE (1,2)) 392 * printf ("This is version 1.2 or later\n") ; 393 * #else 394 * printf ("This is an early version\n") ; 395 * #endif 396 * 397 * Versions 1.1 and earlier of AMD do not include a #define'd version number. 398 */ 399 400 #define AMD_DATE "Jun 20, 2012" 401 #define AMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub)) 402 #define AMD_MAIN_VERSION 2 403 #define AMD_SUB_VERSION 3 404 #define AMD_SUBSUB_VERSION 1 405 #define AMD_VERSION AMD_VERSION_CODE(AMD_MAIN_VERSION,AMD_SUB_VERSION) 406 407 #ifdef __cplusplus 408 } 409 #endif 410 411 #endif 412