old_colamd.c

LU分解求解矩阵方程组的解

    int n_row,
    int n_col,
    int Alen,
    RowInfo Row [],
    ColInfo Col [],
    int A [],
    int head [],
    int n_col2,
    int max_deg,
    int pfree
) ;

PRIVATE void order_children
(
    int n_col,
    ColInfo Col [],
    int p []
) ;

PRIVATE void detect_super_cols
(
#ifndef NDEBUG
    int n_col,
    RowInfo Row [],
#endif
    ColInfo Col [],
    int A [],
    int head [],
    int row_start,
    int row_length
) ;

PRIVATE int garbage_collection
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo Col [],
    int A [],
    int *pfree
) ;

PRIVATE int clear_mark
(
    int n_row,
    RowInfo Row []
) ;

/* ========================================================================== */
/* === Debugging definitions ================================================ */
/* ========================================================================== */

#ifndef NDEBUG

/* === With debugging ======================================================= */

/* stdlib.h: for getenv and atoi, to get debugging level from environment */
#include <stdlib.h>

/* stdio.h:  for printf (no printing if debugging is turned off) */
#include <stdio.h>

PRIVATE void debug_deg_lists
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo Col [],
    int head [],
    int min_score,
    int should,
    int max_deg
) ;

PRIVATE void debug_mark
(
    int n_row,
    RowInfo Row [],
    int tag_mark,
    int max_mark
) ;

PRIVATE void debug_matrix
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo Col [],
    int A []
) ;

PRIVATE void debug_structures
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo Col [],
    int A [],
    int n_col2
) ;

/* the following is the *ONLY* global variable in this file, and is only */
/* present when debugging */

PRIVATE int debug_colamd ;	/* debug print level */

#define DEBUG0(params) { (void) printf params ; }
#define DEBUG1(params) { if (debug_colamd >= 1) (void) printf params ; }
#define DEBUG2(params) { if (debug_colamd >= 2) (void) printf params ; }
#define DEBUG3(params) { if (debug_colamd >= 3) (void) printf params ; }
#define DEBUG4(params) { if (debug_colamd >= 4) (void) printf params ; }

#else

/* === No debugging ========================================================= */

#define DEBUG0(params) ;
#define DEBUG1(params) ;
#define DEBUG2(params) ;
#define DEBUG3(params) ;
#define DEBUG4(params) ;

#endif

/* ========================================================================== */


/* ========================================================================== */
/* === USER-CALLABLE ROUTINES: ============================================== */
/* ========================================================================== */


/* ========================================================================== */
/* === colamd_recommended =================================================== */
/* ========================================================================== */

/*
    The colamd_recommended routine returns the suggested size for Alen.  This
    value has been determined to provide good balance between the number of
    garbage collections and the memory requirements for colamd.
*/

PUBLIC int colamd_recommended	/* returns recommended value of Alen. */
(
    /* === Parameters ======================================================= */

    int nnz,			/* number of nonzeros in A */
    int n_row,			/* number of rows in A */
    int n_col			/* number of columns in A */
)
{
    /* === Local variables ================================================== */

    int minimum ;		/* bare minimum requirements */
    int recommended ;		/* recommended value of Alen */

    if (nnz < 0 || n_row < 0 || n_col < 0)
    {
	/* return -1 if any input argument is corrupted */
	DEBUG0 (("colamd_recommended error!")) ;
	DEBUG0 ((" nnz: %d, n_row: %d, n_col: %d\n", nnz, n_row, n_col)) ;
	return (-1) ;
    }

    minimum =
	2 * (nnz)		/* for A */
	+ (((n_col) + 1) * sizeof (ColInfo) / sizeof (int))	/* for Col */
	+ (((n_row) + 1) * sizeof (RowInfo) / sizeof (int))	/* for Row */
	+ n_col			/* minimum elbow room to guarrantee success */
	+ COLAMD_STATS ;	/* for output statistics */

    /* recommended is equal to the minumum plus enough memory to keep the */
    /* number garbage collections low */
    recommended = minimum + nnz/5 ;

    return (recommended) ;
}


/* ========================================================================== */
/* === colamd_set_defaults ================================================== */
/* ========================================================================== */

/*
    The colamd_set_defaults routine sets the default values of the user-
    controllable parameters for colamd:

	knobs [0]	rows with knobs[0]*n_col entries or more are removed
			prior to ordering.

	knobs [1]	columns with knobs[1]*n_row entries or more are removed
			prior to ordering, and placed last in the column
			permutation.

	knobs [2..19]	unused, but future versions might use this
*/

PUBLIC void colamd_set_defaults
(
    /* === Parameters ======================================================= */

    double knobs [COLAMD_KNOBS]		/* knob array */
)
{
    /* === Local variables ================================================== */

    int i ;

    if (!knobs)
    {
	return ;			/* no knobs to initialize */
    }
    for (i = 0 ; i < COLAMD_KNOBS ; i++)
    {
	knobs [i] = 0 ;
    }
    knobs [COLAMD_DENSE_ROW] = 0.5 ;	/* ignore rows over 50% dense */
    knobs [COLAMD_DENSE_COL] = 0.5 ;	/* ignore columns over 50% dense */
}


/* ========================================================================== */
/* === colamd =============================================================== */
/* ========================================================================== */

/*
    The colamd routine computes a column ordering Q of a sparse matrix
    A such that the LU factorization P(AQ) = LU remains sparse, where P is
    selected via partial pivoting.   The routine can also be viewed as
    providing a permutation Q such that the Cholesky factorization
    (AQ)'(AQ) = LL' remains sparse.

    On input, the nonzero patterns of the columns of A are stored in the
    array A, in order 0 to n_col-1.  A is held in 0-based form (rows in the
    range 0 to n_row-1 and columns in the range 0 to n_col-1).  Row indices
    for column c are located in A [(p [c]) ... (p [c+1]-1)], where p [0] = 0,
    and thus p [n_col] is the number of entries in A.  The matrix is
    destroyed on output.  The row indices within each column do not have to
    be sorted (from small to large row indices), and duplicate row indices
    may be present.  However, colamd will work a little faster if columns are
    sorted and no duplicates are present.  Matlab 5.2 always passes the matrix
    with sorted columns, and no duplicates.

    The integer array A is of size Alen.  Alen must be at least of size
    (where nnz is the number of entries in A):

	nnz			for the input column form of A
	+ nnz			for a row form of A that colamd generates
	+ 6*(n_col+1)		for a ColInfo Col [0..n_col] array
				(this assumes sizeof (ColInfo) is 6 int's).
	+ 4*(n_row+1)		for a RowInfo Row [0..n_row] array
				(this assumes sizeof (RowInfo) is 4 int's).
	+ elbow_room		must be at least n_col.  We recommend at least
				nnz/5 in addition to that.  If sufficient,
				changes in the elbow room affect the ordering
				time only, not the ordering itself.
	+ COLAMD_STATS		for the output statistics

    Colamd returns FALSE is memory is insufficient, or TRUE otherwise.

    On input, the caller must specify:

	n_row			the number of rows of A
	n_col			the number of columns of A
	Alen			the size of the array A
	A [0 ... nnz-1]		the row indices, where nnz = p [n_col]
	A [nnz ... Alen-1]	(need not be initialized by the user)
	p [0 ... n_col]		the column pointers,  p [0] = 0, and p [n_col]
				is the number of entries in A.  Column c of A
				is stored in A [p [c] ... p [c+1]-1].
	knobs [0 ... 19]	a set of parameters that control the behavior
				of colamd.  If knobs is a NULL pointer the
				defaults are used.  The user-callable
				colamd_set_defaults routine sets the default
				parameters.  See that routine for a description
				of the user-controllable parameters.

    If the return value of Colamd is TRUE, then on output:

	p [0 ... n_col-1]	the column permutation. p [0] is the first
				column index, and p [n_col-1] is the last.
				That is, p [k] = j means that column j of A
				is the kth column of AQ.

	A			is undefined on output (the matrix pattern is
				destroyed), except for the following statistics:

	A [0]			the number of dense (or empty) rows ignored
	A [1]			the number of dense (or empty) columms.  These
				are ordered last, in their natural order.
	A [2]			the number of garbage collections performed.
				If this is excessive, then you would have
				gotten your results faster if Alen was larger.
	A [3]			0, if all row indices in each column were in
				sorted order and no duplicates were present.
				1, if there were unsorted or duplicate row
				indices in the input.  You would have gotten
				your results faster if A [3] was returned as 0.

    If the return value of Colamd is FALSE, then A and p are undefined on
    output.
*/

PUBLIC int colamd		/* returns TRUE if successful */
(
    /* === Parameters ======================================================= */

    int n_row,			/* number of rows in A */
    int n_col,			/* number of columns in A */
    int Alen,			/* length of A */
    int A [],			/* row indices of A */
    int p [],			/* pointers to columns in A */
    double knobs [COLAMD_KNOBS]	/* parameters (uses defaults if NULL) */
)
{
    /* === Local variables ================================================== */

    int i ;			/* loop index */
    int nnz ;			/* nonzeros in A */
    int Row_size ;		/* size of Row [], in integers */
    int Col_size ;		/* size of Col [], in integers */
    int elbow_room ;		/* remaining free space */
    RowInfo *Row ;		/* pointer into A of Row [0..n_row] array */
    ColInfo *Col ;		/* pointer into A of Col [0..n_col] array */
    int n_col2 ;		/* number of non-dense, non-empty columns */
    int n_row2 ;		/* number of non-dense, non-empty rows */
    int ngarbage ;		/* number of garbage collections performed */
    int max_deg ;		/* maximum row degree */
    double default_knobs [COLAMD_KNOBS] ;	/* default knobs knobs array */
    int init_result ;		/* return code from initialization */

#ifndef NDEBUG
    debug_colamd = 0 ;		/* no debug printing */
    /* get "D" environment variable, which gives the debug printing level */
    if (getenv ("D")) debug_colamd = atoi (getenv ("D")) ;
    DEBUG0 (("debug version, D = %d (THIS WILL BE SLOOOOW!)\n", debug_colamd)) ;
#endif

    /* === Check the input arguments ======================================== */

    if (n_row < 0 || n_col < 0 || !A || !p)
    {
	/* n_row and n_col must be non-negative, A and p must be present */
	DEBUG0 (("colamd error! %d %d %d\n", n_row, n_col, Alen)) ;
	return (FALSE) ;
    }
    nnz = p [n_col] ;
    if (nnz < 0 || p [0] != 0)
    {
	/* nnz must be non-negative, and p [0] must be zero */
	DEBUG0 (("colamd error! %d %d\n", nnz, p [0])) ;
	return (FALSE) ;
    }

    /* === If no knobs, set default parameters ============================== */

    if (!knobs)
    {
	knobs = default_knobs ;
	colamd_set_defaults (knobs) ;
    }

    /* === Allocate the Row and Col arrays from array A ===================== */

    Col_size = (n_col + 1) * sizeof (ColInfo) / sizeof (int) ;
    Row_size = (n_row + 1) * sizeof (RowInfo) / sizeof (int) ;
    elbow_room = Alen - (2*nnz + Col_size + Row_size) ;
    if (elbow_room < n_col + COLAMD_STATS)
    {
	/* not enough space in array A to perform the ordering */
	DEBUG0 (("colamd error! elbow_room %d, %d\n", elbow_room,n_col)) ;
	return (FALSE) ;
    }
    Alen = 2*nnz + elbow_room ;
    Col  = (ColInfo *) &A [Alen] ;
    Row  = (RowInfo *) &A [Alen + Col_size] ;

    /* === Construct the row and column data structures ===================== */

    init_result = init_rows_cols (n_row, n_col, Row, Col, A, p) ;
    if (init_result == -1)
    {
	/* input matrix is invalid */
	DEBUG0 (("colamd error! matrix invalid\n")) ;
	return (FALSE) ;
    }

    /* === Initialize scores, kill dense rows/columns ======================= */

    init_scoring (n_row, n_col, Row, Col, A, p, knobs,
	&n_row2, &n_col2, &max_deg) ;

    /* === Order the supercolumns =========================================== */

    ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
	n_col2, max_deg, 2*nnz) ;

    /* === Order the non-principal columns ================================== */

    order_children (n_col, Col, p) ;

    /* === Return statistics in A =========================================== */

    for (i = 0 ; i < COLAMD_STATS ; i++)