amd.h 18 KB

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