colamd.c

LU矩阵分解单机版最新版本

		On output, if symamd returns TRUE, the array perm holds the
		permutation P, where perm [0] is the first index in the new
		ordering, and perm [n-1] is the last.  That is, perm [k] = j
		means that row and column j of A is the kth column in PAP',
		where k is in the range 0 to n-1 (perm [0] = j means
		that row and column j of A are the first row and column in
		PAP').  The array is used as a workspace during the ordering,
		which is why it must be of length n+1, not just n.

	    double knobs [COLAMD_KNOBS] ;	Input argument.

		See colamd_set_defaults for a description.

	    int stats [COLAMD_STATS] ;		Output argument.

		Statistics on the ordering, and error status.
		See colamd.h for related definitions.
		Symamd returns FALSE if stats is not present.

		stats [0]:  number of dense or empty row and columns ignored
				(and ordered last in the output permutation 
				perm).  Note that a row/column can become
				"empty" if it contains only "dense" and/or
				"empty" columns/rows.

		stats [1]:  (same as stats [0])

		stats [2]:  number of garbage collections performed.

		stats [3]:  status code.  < 0 is an error code.
			    > 1 is a warning or notice.

			0	OK.  Each column of the input matrix contained
				row indices in increasing order, with no
				duplicates.

			1	OK, but columns of input matrix were jumbled
				(unsorted columns or duplicate entries).  Symamd
				had to do some extra work to sort the matrix
				first and remove duplicate entries, but it
				still was able to return a valid permutation
				(return value of symamd was TRUE).

					stats [4]: highest numbered column that
						is unsorted or has duplicate
						entries.
					stats [5]: last seen duplicate or
						unsorted row index.
					stats [6]: number of duplicate or
						unsorted row indices.

			-1	A is a null pointer

			-2	p is a null pointer

			-3	(unused, see colamd.c)

			-4 	n is negative

					stats [4]: n

			-5	number of nonzeros in matrix is negative

					stats [4]: # of nonzeros (p [n]).

			-6	p [0] is nonzero

					stats [4]: p [0]

			-7	(unused)

			-8	a column has a negative number of entries

					stats [4]: column with < 0 entries
					stats [5]: number of entries in col

			-9	a row index is out of bounds

					stats [4]: column with bad row index
					stats [5]: bad row index
					stats [6]: n_row, # of rows of matrx

			-10	out of memory (unable to allocate temporary
				workspace for M or count arrays using the
				"allocate" routine passed into symamd).

			-999	internal error.  colamd failed to order the
				matrix M, when it should have succeeded.  This
				indicates a bug.  If this (and *only* this)
				error code occurs, please contact the authors.
				Don't contact the authors if you get any other
				error code.

		Future versions may return more statistics in the stats array.

	    void * (*allocate) (size_t, size_t)

	    	A pointer to a function providing memory allocation.  The
		allocated memory must be returned initialized to zero.  For a
		C application, this argument should normally be a pointer to
		calloc.  For a MATLAB mexFunction, the routine mxCalloc is
		passed instead.

	    void (*release) (size_t, size_t)

	    	A pointer to a function that frees memory allocated by the
		memory allocation routine above.  For a C application, this
		argument should normally be a pointer to free.  For a MATLAB
		mexFunction, the routine mxFree is passed instead.


    ----------------------------------------------------------------------------
    colamd_report:
    ----------------------------------------------------------------------------

	C syntax:

	    #include "colamd.h"
	    colamd_report (int stats [COLAMD_STATS]) ;

	Purpose:

	    Prints the error status and statistics recorded in the stats
	    array on the standard error output (for a standard C routine)
	    or on the MATLAB output (for a mexFunction).

	Arguments:

	    int stats [COLAMD_STATS] ;	Input only.  Statistics from colamd.


    ----------------------------------------------------------------------------
    symamd_report:
    ----------------------------------------------------------------------------

	C syntax:

	    #include "colamd.h"
	    symamd_report (int stats [COLAMD_STATS]) ;

	Purpose:

	    Prints the error status and statistics recorded in the stats
	    array on the standard error output (for a standard C routine)
	    or on the MATLAB output (for a mexFunction).

	Arguments:

	    int stats [COLAMD_STATS] ;	Input only.  Statistics from symamd.


*/

/* ========================================================================== */
/* === Scaffolding code definitions  ======================================== */
/* ========================================================================== */

/* Ensure that debugging is turned off: */
#ifndef NDEBUG
#define NDEBUG
#endif /* NDEBUG */

/*
   Our "scaffolding code" philosophy:  In our opinion, well-written library
   code should keep its "debugging" code, and just normally have it turned off
   by the compiler so as not to interfere with performance.  This serves
   several purposes:

   (1) assertions act as comments to the reader, telling you what the code
	expects at that point.  All assertions will always be true (unless
	there really is a bug, of course).

   (2) leaving in the scaffolding code assists anyone who would like to modify
	the code, or understand the algorithm (by reading the debugging output,
	one can get a glimpse into what the code is doing).

   (3) (gasp!) for actually finding bugs.  This code has been heavily tested
	and "should" be fully functional and bug-free ... but you never know...

    To enable debugging, comment out the "#define NDEBUG" above.  For a MATLAB
    mexFunction, you will also need to modify mexopts.sh to remove the -DNDEBUG
    definition.  The code will become outrageously slow when debugging is
    enabled.  To control the level of debugging output, set an environment
    variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
    you should see the following message on the standard output:

    	colamd: debug version, D = 1 (THIS WILL BE SLOW!)

    or a similar message for symamd.  If you don't, then debugging has not
    been enabled.

*/

/* ========================================================================== */
/* === Include files ======================================================== */
/* ========================================================================== */

#include "colamd.h"
#include <limits.h>

#ifdef MATLAB_MEX_FILE
#include "mex.h"
#include "matrix.h"
#else
#include <stdio.h>
#include <assert.h>
#endif /* MATLAB_MEX_FILE */

/* ========================================================================== */
/* === Definitions ========================================================== */
/* ========================================================================== */

/* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
#define PUBLIC
#define PRIVATE static

#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define MIN(a,b) (((a) < (b)) ? (a) : (b))

#define ONES_COMPLEMENT(r) (-(r)-1)

/* -------------------------------------------------------------------------- */
/* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
/* -------------------------------------------------------------------------- */

#ifndef TRUE
#define TRUE (1)
#endif

#ifndef FALSE
#define FALSE (0)
#endif

/* -------------------------------------------------------------------------- */

#define EMPTY	(-1)

/* Row and column status */
#define ALIVE	(0)
#define DEAD	(-1)

/* Column status */
#define DEAD_PRINCIPAL		(-1)
#define DEAD_NON_PRINCIPAL	(-2)

/* Macros for row and column status update and checking. */
#define ROW_IS_DEAD(r)			ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
#define ROW_IS_MARKED_DEAD(row_mark)	(row_mark < ALIVE)
#define ROW_IS_ALIVE(r)			(Row [r].shared2.mark >= ALIVE)
#define COL_IS_DEAD(c)			(Col [c].start < ALIVE)
#define COL_IS_ALIVE(c)			(Col [c].start >= ALIVE)
#define COL_IS_DEAD_PRINCIPAL(c)	(Col [c].start == DEAD_PRINCIPAL)
#define KILL_ROW(r)			{ Row [r].shared2.mark = DEAD ; }
#define KILL_PRINCIPAL_COL(c)		{ Col [c].start = DEAD_PRINCIPAL ; }
#define KILL_NON_PRINCIPAL_COL(c)	{ Col [c].start = DEAD_NON_PRINCIPAL ; }

/* ========================================================================== */
/* === Colamd reporting mechanism =========================================== */
/* ========================================================================== */

#ifdef MATLAB_MEX_FILE

/* use mexPrintf in a MATLAB mexFunction, for debugging and statistics output */
#define PRINTF mexPrintf

/* In MATLAB, matrices are 1-based to the user, but 0-based internally */
#define INDEX(i) ((i)+1)

#else

/* Use printf in standard C environment, for debugging and statistics output. */
/* Output is generated only if debugging is enabled at compile time, or if */
/* the caller explicitly calls colamd_report or symamd_report. */
#define PRINTF printf

/* In C, matrices are 0-based and indices are reported as such in *_report */
#define INDEX(i) (i)

#endif /* MATLAB_MEX_FILE */

/* ========================================================================== */
/* === Prototypes of PRIVATE routines ======================================= */
/* ========================================================================== */

PRIVATE int init_rows_cols
(
    int n_row,
    int n_col,
    Colamd_Row Row [],
    Colamd_Col Col [],
    int A [],
    int p [],
    int stats [COLAMD_STATS]
) ;

PRIVATE void init_scoring
(
    int n_row,
    int n_col,
    Colamd_Row Row [],
    Colamd_Col Col [],
    int A [],
    int head [],
    double knobs [COLAMD_KNOBS],
    int *p_n_row2,
    int *p_n_col2,
    int *p_max_deg
) ;

PRIVATE int find_ordering
(
    int n_row,
    int n_col,
    int Alen,
    Colamd_Row Row [],
    Colamd_Col Col [],
    int A [],
    int head [],
    int n_col2,
    int max_deg,
    int pfree
) ;

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

PRIVATE void detect_super_cols
(

#ifndef NDEBUG
    int n_col,
    Colamd_Row Row [],
#endif /* NDEBUG */

    Colamd_Col Col [],
    int A [],
    int head [],
    int row_start,
    int row_length
) ;

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

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

PRIVATE void print_report
(
    char *method,
    int stats [COLAMD_STATS]
) ;

/* ========================================================================== */
/* === Debugging prototypes and definitions ================================= */
/* ========================================================================== */

#ifndef NDEBUG

/* colamd_debug is the *ONLY* global variable, and is only */
/* present when debugging */

PRIVATE int colamd_debug ;	/* debug print level */

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

#ifdef MATLAB_MEX_FILE
#define ASSERT(expression) (mxAssert ((expression), ""))
#else
#define ASSERT(expression) (assert (expression))
#endif /* MATLAB_MEX_FILE */

PRIVATE void colamd_get_debug	/* gets the debug print level from getenv */
(
    char *method
) ;

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

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

PRIVATE void debug_matrix
(
    int n_row,
    int n_col,
    Colamd_Row Row [],
    Colamd_Col Col [],
    int A []
) ;

PRIVATE void debug_structures
(