old_colamd.c

SuperLU 2.2版本。对大型、稀疏、非对称的线性系统的直接求解

/* ========================================================================== */
/* === colamd - a sparse matrix column ordering algorithm =================== */
/* ========================================================================== */

/*
    colamd:  An approximate minimum degree column ordering algorithm.

    Purpose:

	Colamd computes a permutation Q such that the Cholesky factorization of
	(AQ)'(AQ) has less fill-in and requires fewer floating point operations
	than A'A.  This also provides a good ordering for sparse partial
	pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
	factorization, and P is computed during numerical factorization via
	conventional partial pivoting with row interchanges.  Colamd is the
	column ordering method used in SuperLU, part of the ScaLAPACK library.
	It is also available as user-contributed software for Matlab 5.2,
	available from MathWorks, Inc. (http://www.mathworks.com).  This
	routine can be used in place of COLMMD in Matlab.  By default, the \
	and / operators in Matlab perform a column ordering (using COLMMD)
	prior to LU factorization using sparse partial pivoting, in the
	built-in Matlab LU(A) routine.

    Authors:

	The authors of the code itself are Stefan I. Larimore and Timothy A.
	Davis ([e-mail]), University of Florida.  The algorithm was
	developed in collaboration with John Gilbert, Xerox PARC, and Esmond
	Ng, Oak Ridge National Laboratory.

    Date:

	August 3, 1998.  Version 1.0.

    Acknowledgements:

	This work was supported by the National Science Foundation, under
	grants DMS-9504974 and DMS-9803599.

    Notice:

	Copyright (c) 1998 by the University of Florida.  All Rights Reserved.

	THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
	EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.

	Permission is hereby granted to use or copy this program for any
	purpose, provided the above notices are retained on all copies.
	User documentation of any code that uses this code must cite the
	Authors, the Copyright, and "Used by permission."  If this code is
	accessible from within Matlab, then typing "help colamd" or "colamd"
	(with no arguments) must cite the Authors.  Permission to modify the
	code and to distribute modified code is granted, provided the above
	notices are retained, and a notice that the code was modified is
	included with the above copyright notice.  You must also retain the
	Availability information below, of the original version.

	This software is provided free of charge.

    Availability:

	This file is located at

		http://www.cise.ufl.edu/~davis/colamd/colamd.c

	The colamd.h file is required, located in the same directory.
	The colamdmex.c file provides a Matlab interface for colamd.
	The symamdmex.c file provides a Matlab interface for symamd, which is
	a symmetric ordering based on this code, colamd.c.  All codes are
	purely ANSI C compliant (they use no Unix-specific routines, include
	files, etc.).
*/

/* ========================================================================== */
/* === Description of user-callable routines ================================ */
/* ========================================================================== */

/*
    Each user-callable routine (declared as PUBLIC) is briefly described below.
    Refer to the comments preceding each routine for more details.

    ----------------------------------------------------------------------------
    colamd_recommended:
    ----------------------------------------------------------------------------

	Usage:

	    Alen = colamd_recommended (nnz, n_row, n_col) ;

	Purpose:

	    Returns recommended value of Alen for use by colamd.  Returns -1
	    if any input argument is negative.

	Arguments:

	    int nnz ;		Number of nonzeros in the matrix A.  This must
				be the same value as p [n_col] in the call to
				colamd - otherwise you will get a wrong value
				of the recommended memory to use.
	    int n_row ;		Number of rows in the matrix A.
	    int n_col ;		Number of columns in the matrix A.

    ----------------------------------------------------------------------------
    colamd_set_defaults:
    ----------------------------------------------------------------------------

	Usage:

	    colamd_set_defaults (knobs) ;

	Purpose:

	    Sets the default parameters.

	Arguments:

	    double knobs [COLAMD_KNOBS] ;	Output only.

		Rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) entries
		are removed prior to ordering.  Columns with more than
		(knobs [COLAMD_DENSE_COL] * n_row) entries are removed
		prior to ordering, and placed last in the output column
		ordering.  Default values of these two knobs are both 0.5.
		Currently, only knobs [0] and knobs [1] are used, but future
		versions may use more knobs.  If so, they will be properly set
		to their defaults by the future version of colamd_set_defaults,
		so that the code that calls colamd will not need to change,
		assuming that you either use colamd_set_defaults, or pass a
		(double *) NULL pointer as the knobs array to colamd.

    ----------------------------------------------------------------------------
    colamd:
    ----------------------------------------------------------------------------

	Usage:

	    colamd (n_row, n_col, Alen, A, p, knobs) ;

	Purpose:

	    Computes a column ordering (Q) of A such that P(AQ)=LU or
	    (AQ)'AQ=LL' have less fill-in and require fewer floating point
	    operations than factorizing the unpermuted matrix A or A'A,
	    respectively.

	Arguments:

	    int n_row ;

		Number of rows in the matrix A.
		Restriction:  n_row >= 0.
		Colamd returns FALSE if n_row is negative.

	    int n_col ;

		Number of columns in the matrix A.
		Restriction:  n_col >= 0.
		Colamd returns FALSE if n_col is negative.

	    int Alen ;

		Restriction (see note):
		Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col + COLAMD_STATS
		Colamd returns FALSE if these conditions are not met.

		Note:  this restriction makes an modest assumption regarding
		the size of the two typedef'd structures, below.  We do,
		however, guarantee that
		Alen >= colamd_recommended (nnz, n_row, n_col)
		will be sufficient.

	    int A [Alen] ;	Input argument, stats on output.

		A is an integer array of size Alen.  Alen must be at least as
		large as the bare minimum value given above, but this is very
		low, and can result in excessive run time.  For best
		performance, we recommend that Alen be greater than or equal to
		colamd_recommended (nnz, n_row, n_col), which adds
		nnz/5 to the bare minimum value given above.

		On input, the row indices of the entries in column c of the
		matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
		in a given column c need not be in ascending order, and
		duplicate row indices may be be present.  However, colamd will
		work a little faster if both of these conditions are met
		(Colamd puts the matrix into this format, if it finds that the
		the conditions are not met).

		The matrix is 0-based.  That is, rows are in the range 0 to
		n_row-1, and columns are in the range 0 to n_col-1.  Colamd
		returns FALSE if any row index is out of range.

		The contents of A are modified during ordering, and are thus
		undefined on output with the exception of a few statistics
		about the ordering (A [0..COLAMD_STATS-1]):
		A [0]:  number of dense or empty rows ignored.
		A [1]:  number of dense or empty columns ignored (and ordered
			last in the output permutation p)
		A [2]:  number of garbage collections performed.
		A [3]:  0, if all row indices in each column were in sorted
			  order, and no duplicates were present.
			1, otherwise (in which case colamd had to do more work)
		Note that a row can become "empty" if it contains only
		"dense" and/or "empty" columns, and similarly a column can
		become "empty" if it only contains "dense" and/or "empty" rows.
		Future versions may return more statistics in A, but the usage
		of these 4 entries in A will remain unchanged.

	    int p [n_col+1] ;	Both input and output argument.

		p is an integer array of size n_col+1.  On input, it holds the
		"pointers" for the column form of the matrix A.  Column c of
		the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
		entry, p [0], must be zero, and p [c] <= p [c+1] must hold
		for all c in the range 0 to n_col-1.  The value p [n_col] is
		thus the total number of entries in the pattern of the matrix A.
		Colamd returns FALSE if these conditions are not met.

		On output, if colamd returns TRUE, the array p holds the column
		permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
		the first column index in the new ordering, and p [n_col-1] is
		the last.  That is, p [k] = j means that column j of A is the
		kth pivot column, in AQ, where k is in the range 0 to n_col-1
		(p [0] = j means that column j of A is the first column in AQ).

		If colamd returns FALSE, then no permutation is returned, and
		p is undefined on output.

	    double knobs [COLAMD_KNOBS] ;	Input only.

		See colamd_set_defaults for a description.  If the knobs array
		is not present (that is, if a (double *) NULL pointer is passed
		in its place), then the default values of the parameters are
		used instead.

*/


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

/* limits.h:  the largest positive integer (INT_MAX) */
#include <limits.h>

/* colamd.h:  knob array size, stats output size, and global prototypes */
#include "colamd.h"

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

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

/* assert.h:  the assert macro (no debugging if NDEBUG is defined) */
#include <assert.h>

/*
   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.  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).
*/

/* ========================================================================== */
/* === Row and Column structures ============================================ */
/* ========================================================================== */

typedef struct ColInfo_struct
{
    int start ;		/* index for A of first row in this column, or DEAD */
			/* if column is dead */
    int length ;	/* number of rows in this column */
    union
    {
	int thickness ;	/* number of original columns represented by this */
			/* col, if the column is alive */
	int parent ;	/* parent in parent tree super-column structure, if */
			/* the column is dead */
    } shared1 ;
    union
    {
	int score ;	/* the score used to maintain heap, if col is alive */
	int order ;	/* pivot ordering of this column, if col is dead */
    } shared2 ;
    union
    {
	int headhash ;	/* head of a hash bucket, if col is at the head of */
			/* a degree list */
	int hash ;	/* hash value, if col is not in a degree list */
	int prev ;	/* previous column in degree list, if col is in a */
			/* degree list (but not at the head of a degree list) */
    } shared3 ;
    union
    {
	int degree_next ;	/* next column, if col is in a degree list */
	int hash_next ;		/* next column, if col is in a hash list */
    } shared4 ;

} ColInfo ;

typedef struct RowInfo_struct
{
    int start ;		/* index for A of first col in this row */
    int length ;	/* number of principal columns in this row */
    union
    {
	int degree ;	/* number of principal & non-principal columns in row */
	int p ;		/* used as a row pointer in init_rows_cols () */
    } shared1 ;
    union
    {
	int mark ;	/* for computing set differences and marking dead rows*/
	int first_column ;/* first column in row (used in garbage collection) */
    } shared2 ;

} RowInfo ;

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

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

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

#define TRUE	(1)
#define FALSE	(0)
#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 ; }

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

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

PRIVATE int init_rows_cols
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo Col [],
    int A [],
    int p []
) ;

PRIVATE void init_scoring
(
    int n_row,
    int n_col,
    RowInfo Row [],
    ColInfo 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
(