amd.h

数学算法的实现库。可以实现常见的线性计算。

/* ========================================================================= */
/* === AMD:  approximate minimum degree ordering =========================== */
/* ========================================================================= */

/* ------------------------------------------------------------------------- */
/* AMD Version 1.1 (Jan. 21, 2004), Copyright (c) 2004 by Timothy A. Davis,  */
/* Patrick R. Amestoy, and Iain S. Duff.  See ../README for License.         */
/* email: [e-mail]    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 "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

int amd_order (		    /* returns 0 if OK, negative value if error */
    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 */
) ;

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

/* Input arguments (not modified):
 *
 *	n: the matrix A is n-by-n.
 *	Ap: an int/long array of size n+1, containing the column pointers of A.
 *	Ai: an int/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/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].
 * The row indices must appear in ascending order in each column, and there
 * must not be any duplicate entries.  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, or use the amd_preprocess
 * routine discussed below.
 *
 * 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.  For all j in the range 0 to n-1,
 *	and for all p in the range Ap [j] to Ap [j+1]-2, Ai [p] < Ai [p+1] must
 *	hold.  Ai [0..nz-1] must be in the range 0 to n-1.  To avoid integer
 *	overflow, (2.4*nz + 8*n) < INT_MAX / sizeof (int) for must hold for the
 *	"int" version. (2.4*nz + 8*n) < LONG_MAX / sizeof (long) must hold
 *	for the "long" version.  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.
 *
 * 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_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.
 *