amd.h 16.9 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
/* ========================================================================= */
/* === AMD:  approximate minimum degree ordering =========================== */
/* ========================================================================= */

/* ------------------------------------------------------------------------- */
/* AMD Version 2.2, Copyright (c) 2007 by Timothy A. Davis,                  */
/* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
/* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
/* web: http://www.cise.ufl.edu/research/sparse/amd                          */
/* ------------------------------------------------------------------------- */

/* AMD finds a symmetric ordering P of a matrix A so that the Cholesky
 * factorization of P*A*P' has fewer nonzeros and takes less work than the
 * Cholesky factorization of A.  If A is not symmetric, then it performs its
 * ordering on the matrix A+A'.  Two sets of user-callable routines are
 * provided, one for int integers and the other for UF_long integers.
 *
 * The method is based on the approximate minimum degree algorithm, discussed
 * in Amestoy, Davis, and Duff, "An approximate degree ordering algorithm",
 * SIAM Journal of Matrix Analysis and Applications, vol. 17, no. 4, pp.
 * 886-905, 1996.  This package can perform both the AMD ordering (with
 * aggressive absorption), and the AMDBAR ordering (without aggressive
 * absorption) discussed in the above paper.  This package differs from the
 * Fortran codes discussed in the paper:
 *
 *	(1) it can ignore "dense" rows and columns, leading to faster run times
 *	(2) it computes the ordering of A+A' if A is not symmetric
 *	(3) it is followed by a depth-first post-ordering of the assembly tree
 *	    (or supernodal elimination tree)
 *
 * For historical reasons, the Fortran versions, amd.f and amdbar.f, have
 * been left (nearly) unchanged.  They compute the identical ordering as
 * described in the above paper.
 */

#ifndef AMD_H
#define AMD_H

/* make it easy for C++ programs to include AMD */
#ifdef __cplusplus
extern "C" {
#endif

/* get the definition of size_t: */
#include <stddef.h>

/* define UF_long */
#include "UFconfig.h"

int amd_order		    /* returns AMD_OK, AMD_OK_BUT_JUMBLED,
			     * AMD_INVALID, or AMD_OUT_OF_MEMORY */
(
    int n,		    /* A is n-by-n.  n must be >= 0. */
    const int Ap [ ],	    /* column pointers for A, of size n+1 */
    const int Ai [ ],	    /* row indices of A, of size nz = Ap [n] */
    int P [ ],		    /* output permutation, of size n */
    double Control [ ],	    /* input Control settings, of size AMD_CONTROL */
    double Info [ ]	    /* output Info statistics, of size AMD_INFO */
) ;

UF_long amd_l_order	    /* see above for description of arguments */
(
    UF_long n,
    const UF_long Ap [ ],
    const UF_long Ai [ ],
    UF_long P [ ],
    double Control [ ],
    double Info [ ]
) ;

/* Input arguments (not modified):
 *
 *	n: the matrix A is n-by-n.
 *	Ap: an int/UF_long array of size n+1, containing column pointers of A.
 *	Ai: an int/UF_long array of size nz, containing the row indices of A,
 *	    where nz = Ap [n].
 *	Control:  a double array of size AMD_CONTROL, containing control
 *	    parameters.  Defaults are used if Control is NULL.
 *
 * Output arguments (not defined on input):
 *
 *	P: an int/UF_long array of size n, containing the output permutation. If
 *	    row i is the kth pivot row, then P [k] = i.  In MATLAB notation,
 *	    the reordered matrix is A (P,P).
 *	Info: a double array of size AMD_INFO, containing statistical
 *	    information.  Ignored if Info is NULL.
 *
 * On input, the matrix A is stored in column-oriented form.  The row indices
 * of nonzero entries in column j are stored in Ai [Ap [j] ... Ap [j+1]-1].
 *
 * If the row indices appear in ascending order in each column, and there
 * are no duplicate entries, then amd_order is slightly more efficient in
 * terms of time and memory usage.  If this condition does not hold, a copy
 * of the matrix is created (where these conditions do hold), and the copy is
 * ordered.  This feature is new to v2.0 (v1.2 and earlier required this
 * condition to hold for the input matrix).
 * 
 * Row indices must be in the range 0 to
 * n-1.  Ap [0] must be zero, and thus nz = Ap [n] is the number of nonzeros
 * in A.  The array Ap is of size n+1, and the array Ai is of size nz = Ap [n].
 * The matrix does not need to be symmetric, and the diagonal does not need to
 * be present (if diagonal entries are present, they are ignored except for
 * the output statistic Info [AMD_NZDIAG]).  The arrays Ai and Ap are not
 * modified.  This form of the Ap and Ai arrays to represent the nonzero
 * pattern of the matrix A is the same as that used internally by MATLAB.
 * If you wish to use a more flexible input structure, please see the
 * umfpack_*_triplet_to_col routines in the UMFPACK package, at
 * http://www.cise.ufl.edu/research/sparse/umfpack.
 *
 * Restrictions:  n >= 0.  Ap [0] = 0.  Ap [j] <= Ap [j+1] for all j in the
 *	range 0 to n-1.  nz = Ap [n] >= 0.  Ai [0..nz-1] must be in the range 0
 *	to n-1.  Finally, Ai, Ap, and P must not be NULL.  If any of these
 *	restrictions are not met, AMD returns AMD_INVALID.
 *
 * AMD returns:
 *
 *	AMD_OK if the matrix is valid and sufficient memory can be allocated to
 *	    perform the ordering.
 *
 *	AMD_OUT_OF_MEMORY if not enough memory can be allocated.
 *
 *	AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if P is
 *	    NULL.
 *
 *	AMD_OK_BUT_JUMBLED if the matrix had unsorted columns, and/or duplicate
 *	    entries, but was otherwise valid.
 *
 * The AMD routine first forms the pattern of the matrix A+A', and then
 * computes a fill-reducing ordering, P.  If P [k] = i, then row/column i of
 * the original is the kth pivotal row.  In MATLAB notation, the permuted
 * matrix is A (P,P), except that 0-based indexing is used instead of the
 * 1-based indexing in MATLAB.
 *
 * The Control array is used to set various parameters for AMD.  If a NULL
 * pointer is passed, default values are used.  The Control array is not
 * modified.
 *
 *	Control [AMD_DENSE]:  controls the threshold for "dense" rows/columns.
 *	    A dense row/column in A+A' can cause AMD to spend a lot of time in
 *	    ordering the matrix.  If Control [AMD_DENSE] >= 0, rows/columns
 *	    with more than Control [AMD_DENSE] * sqrt (n) entries are ignored
 *	    during the ordering, and placed last in the output order.  The
 *	    default value of Control [AMD_DENSE] is 10.  If negative, no
 *	    rows/columns are treated as "dense".  Rows/columns with 16 or
 *	    fewer off-diagonal entries are never considered "dense".
 *
 *	Control [AMD_AGGRESSIVE]: controls whether or not to use aggressive
 *	    absorption, in which a prior element is absorbed into the current
 *	    element if is a subset of the current element, even if it is not
 *	    adjacent to the current pivot element (refer to Amestoy, Davis,
 *	    & Duff, 1996, for more details).  The default value is nonzero,
 *	    which means to perform aggressive absorption.  This nearly always
 *	    leads to a better ordering (because the approximate degrees are
 *	    more accurate) and a lower execution time.  There are cases where
 *	    it can lead to a slightly worse ordering, however.  To turn it off,
 *	    set Control [AMD_AGGRESSIVE] to 0.
 *
 *	Control [2..4] are not used in the current version, but may be used in
 *	    future versions.
 *
 * The Info array provides statistics about the ordering on output.  If it is
 * not present, the statistics are not returned.  This is not an error
 * condition.
 * 
 *	Info [AMD_STATUS]:  the return value of AMD, either AMD_OK,
 *	    AMD_OK_BUT_JUMBLED, AMD_OUT_OF_MEMORY, or AMD_INVALID.
 *
 *	Info [AMD_N]: n, the size of the input matrix
 *
 *	Info [AMD_NZ]: the number of nonzeros in A, nz = Ap [n]
 *
 *	Info [AMD_SYMMETRY]:  the symmetry of the matrix A.  It is the number
 *	    of "matched" off-diagonal entries divided by the total number of
 *	    off-diagonal entries.  An entry A(i,j) is matched if A(j,i) is also
 *	    an entry, for any pair (i,j) for which i != j.  In MATLAB notation,
 *		S = spones (A) ;
 *		B = tril (S, -1) + triu (S, 1) ;
 *		symmetry = nnz (B & B') / nnz (B) ;
 *
 *	Info [AMD_NZDIAG]: the number of entries on the diagonal of A.
 *
 *	Info [AMD_NZ_A_PLUS_AT]:  the number of nonzeros in A+A', excluding the
 *	    diagonal.  If A is perfectly symmetric (Info [AMD_SYMMETRY] = 1)
 *	    with a fully nonzero diagonal, then Info [AMD_NZ_A_PLUS_AT] = nz-n
 *	    (the smallest possible value).  If A is perfectly unsymmetric
 *	    (Info [AMD_SYMMETRY] = 0, for an upper triangular matrix, for
 *	    example) with no diagonal, then Info [AMD_NZ_A_PLUS_AT] = 2*nz
 *	    (the largest possible value).
 *
 *	Info [AMD_NDENSE]: the number of "dense" rows/columns of A+A' that were
 *	    removed from A prior to ordering.  These are placed last in the
 *	    output order P.
 *
 *	Info [AMD_MEMORY]: the amount of memory used by AMD, in bytes.  In the
 *	    current version, this is 1.2 * Info  [AMD_NZ_A_PLUS_AT] + 9*n
 *	    times the size of an integer.  This is at most 2.4nz + 9n.  This
 *	    excludes the size of the input arguments Ai, Ap, and P, which have
 *	    a total size of nz + 2*n + 1 integers.
 *
 *	Info [AMD_NCMPA]: the number of garbage collections performed.
 *
 *	Info [AMD_LNZ]: the number of nonzeros in L (excluding the diagonal).
 *	    This is a slight upper bound because mass elimination is combined
 *	    with the approximate degree update.  It is a rough upper bound if
 *	    there are many "dense" rows/columns.  The rest of the statistics,
 *	    below, are also slight or rough upper bounds, for the same reasons.
 *	    The post-ordering of the assembly tree might also not exactly
 *	    correspond to a true elimination tree postordering.
 *
 *	Info [AMD_NDIV]: the number of divide operations for a subsequent LDL'
 *	    or LU factorization of the permuted matrix A (P,P).
 *
 *	Info [AMD_NMULTSUBS_LDL]:  the number of multiply-subtract pairs for a
 *	    subsequent LDL' factorization of A (P,P).
 *
 *	Info [AMD_NMULTSUBS_LU]:  the number of multiply-subtract pairs for a
 *	    subsequent LU factorization of A (P,P), assuming that no numerical
 *	    pivoting is required.
 *
 *	Info [AMD_DMAX]:  the maximum number of nonzeros in any column of L,
 *	    including the diagonal.
 *
 *	Info [14..19] are not used in the current version, but may be used in
 *	    future versions.
 */    

/* ------------------------------------------------------------------------- */
/* direct interface to AMD */
/* ------------------------------------------------------------------------- */

/* amd_2 is the primary AMD ordering routine.  It is not meant to be
 * user-callable because of its restrictive inputs and because it destroys
 * the user's input matrix.  It does not check its inputs for errors, either.
 * However, if you can work with these restrictions it can be faster than
 * amd_order and use less memory (assuming that you can create your own copy
 * of the matrix for AMD to destroy).  Refer to AMD/Source/amd_2.c for a
 * description of each parameter. */

void amd_2
(
    int n,
    int Pe [ ],
    int Iw [ ],
    int Len [ ],
    int iwlen,
    int pfree,
    int Nv [ ],
    int Next [ ], 
    int Last [ ],
    int Head [ ],
    int Elen [ ],
    int Degree [ ],
    int W [ ],
    double Control [ ],
    double Info [ ]
) ;

void amd_l2
(
    UF_long n,
    UF_long Pe [ ],
    UF_long Iw [ ],
    UF_long Len [ ],
    UF_long iwlen,
    UF_long pfree,
    UF_long Nv [ ],
    UF_long Next [ ], 
    UF_long Last [ ],
    UF_long Head [ ],
    UF_long Elen [ ],
    UF_long Degree [ ],
    UF_long W [ ],
    double Control [ ],
    double Info [ ]
) ;

/* ------------------------------------------------------------------------- */
/* amd_valid */
/* ------------------------------------------------------------------------- */

/* Returns AMD_OK or AMD_OK_BUT_JUMBLED if the matrix is valid as input to
 * amd_order; the latter is returned if the matrix has unsorted and/or
 * duplicate row indices in one or more columns.  Returns AMD_INVALID if the
 * matrix cannot be passed to amd_order.  For amd_order, the matrix must also
 * be square.  The first two arguments are the number of rows and the number
 * of columns of the matrix.  For its use in AMD, these must both equal n.
 *
 * NOTE: this routine returned TRUE/FALSE in v1.2 and earlier.
 */

int amd_valid
(
    int n_row,		    /* # of rows */
    int n_col,		    /* # of columns */
    const int Ap [ ],	    /* column pointers, of size n_col+1 */
    const int Ai [ ]	    /* row indices, of size Ap [n_col] */
) ;

UF_long amd_l_valid
(
    UF_long n_row,
    UF_long n_col,
    const UF_long Ap [ ],
    const UF_long Ai [ ]
) ;

/* ------------------------------------------------------------------------- */
/* AMD memory manager and printf routines */
/* ------------------------------------------------------------------------- */

/* The user can redefine these to change the malloc, free, and printf routines
 * that AMD uses. */

#ifndef EXTERN
#define EXTERN extern
#endif

EXTERN void *(*amd_malloc) (size_t) ;		    /* pointer to malloc */
EXTERN void (*amd_free) (void *) ;		    /* pointer to free */
EXTERN void *(*amd_realloc) (void *, size_t) ;	    /* pointer to realloc */
EXTERN void *(*amd_calloc) (size_t, size_t) ;	    /* pointer to calloc */
EXTERN int (*amd_printf) (const char *, ...) ;	    /* pointer to printf */

/* ------------------------------------------------------------------------- */
/* AMD Control and Info arrays */
/* ------------------------------------------------------------------------- */

/* amd_defaults:  sets the default control settings */
void amd_defaults   (double Control [ ]) ;
void amd_l_defaults (double Control [ ]) ;

/* amd_control: prints the control settings */
void amd_control    (double Control [ ]) ;
void amd_l_control  (double Control [ ]) ;

/* amd_info: prints the statistics */
void amd_info       (double Info [ ]) ;
void amd_l_info     (double Info [ ]) ;

#define AMD_CONTROL 5	    /* size of Control array */
#define AMD_INFO 20	    /* size of Info array */

/* contents of Control */
#define AMD_DENSE 0	    /* "dense" if degree > Control [0] * sqrt (n) */
#define AMD_AGGRESSIVE 1    /* do aggressive absorption if Control [1] != 0 */

/* default Control settings */
#define AMD_DEFAULT_DENSE 10.0	    /* default "dense" degree 10*sqrt(n) */
#define AMD_DEFAULT_AGGRESSIVE 1    /* do aggressive absorption by default */

/* contents of Info */
#define AMD_STATUS 0	    /* return value of amd_order and amd_l_order */
#define AMD_N 1		    /* A is n-by-n */
#define AMD_NZ 2	    /* number of nonzeros in A */ 
#define AMD_SYMMETRY 3	    /* symmetry of pattern (1 is sym., 0 is unsym.) */
#define AMD_NZDIAG 4	    /* # of entries on diagonal */
#define AMD_NZ_A_PLUS_AT 5  /* nz in A+A' */
#define AMD_NDENSE 6	    /* number of "dense" rows/columns in A */
#define AMD_MEMORY 7	    /* amount of memory used by AMD */
#define AMD_NCMPA 8	    /* number of garbage collections in AMD */
#define AMD_LNZ 9	    /* approx. nz in L, excluding the diagonal */
#define AMD_NDIV 10	    /* number of fl. point divides for LU and LDL' */
#define AMD_NMULTSUBS_LDL 11 /* number of fl. point (*,-) pairs for LDL' */
#define AMD_NMULTSUBS_LU 12  /* number of fl. point (*,-) pairs for LU */
#define AMD_DMAX 13	     /* max nz. in any column of L, incl. diagonal */

/* ------------------------------------------------------------------------- */
/* return values of AMD */
/* ------------------------------------------------------------------------- */

#define AMD_OK 0		/* success */
#define AMD_OUT_OF_MEMORY -1	/* malloc failed, or problem too large */
#define AMD_INVALID -2		/* input arguments are not valid */
#define AMD_OK_BUT_JUMBLED 1	/* input matrix is OK for amd_order, but
    * columns were not sorted, and/or duplicate entries were present.  AMD had
    * to do extra work before ordering the matrix.  This is a warning, not an
    * error.  */

/* ========================================================================== */
/* === AMD version ========================================================== */
/* ========================================================================== */

/* AMD Version 1.2 and later include the following definitions.
 * As an example, to test if the version you are using is 1.2 or later:
 *
 * #ifdef AMD_VERSION
 *	if (AMD_VERSION >= AMD_VERSION_CODE (1,2)) ...
 * #endif
 *
 * This also works during compile-time:
 *
 *	#if defined(AMD_VERSION) && (AMD_VERSION >= AMD_VERSION_CODE (1,2))
 *	    printf ("This is version 1.2 or later\n") ;
 *	#else
 *	    printf ("This is an early version\n") ;
 *	#endif
 *
 * Versions 1.1 and earlier of AMD do not include a #define'd version number.
 */

#define AMD_DATE "May 31, 2007"
#define AMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define AMD_MAIN_VERSION 2
#define AMD_SUB_VERSION 2
#define AMD_SUBSUB_VERSION 0
#define AMD_VERSION AMD_VERSION_CODE(AMD_MAIN_VERSION,AMD_SUB_VERSION)

#ifdef __cplusplus
}
#endif

#endif