old_colamd.c

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

	    *new_cp++ = row ;
	    /* add row's external degree */
	    score += Row [row].shared1.degree - 1 ;
	    /* guard against integer overflow */
	    score = MIN (score, n_col) ;
	}
	/* determine pruned column length */
	col_length = (int) (new_cp - &A [Col [c].start]) ;
	if (col_length == 0)
	{
	    /* a newly-made null column (all rows in this col are "dense" */
	    /* and have already been killed) */
	    DEBUG0 (("Newly null killed: %d\n", c)) ;
	    Col [c].shared2.order = --n_col2 ;
	    KILL_PRINCIPAL_COL (c) ;
	}
	else
	{
	    /* set column length and set score */
	    assert (score >= 0) ;
	    assert (score <= n_col) ;
	    Col [c].length = col_length ;
	    Col [c].shared2.score = score ;
	}
    }
    DEBUG0 (("Dense, null, and newly-null columns killed: %d\n",n_col-n_col2)) ;

    /* At this point, all empty rows and columns are dead.  All live columns */
    /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
    /* yet).  Rows may contain dead columns, but all live rows contain at */
    /* least one live column. */

#ifndef NDEBUG
    debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
#endif

    /* === Initialize degree lists ========================================== */

#ifndef NDEBUG
    debug_count = 0 ;
#endif

    /* clear the hash buckets */
    for (c = 0 ; c <= n_col ; c++)
    {
	head [c] = EMPTY ;
    }
    min_score = n_col ;
    /* place in reverse order, so low column indices are at the front */
    /* of the lists.  This is to encourage natural tie-breaking */
    for (c = n_col-1 ; c >= 0 ; c--)
    {
	/* only add principal columns to degree lists */
	if (COL_IS_ALIVE (c))
	{
	    DEBUG4 (("place %d score %d minscore %d ncol %d\n",
		c, Col [c].shared2.score, min_score, n_col)) ;

	    /* === Add columns score to DList =============================== */

	    score = Col [c].shared2.score ;

	    assert (min_score >= 0) ;
	    assert (min_score <= n_col) ;
	    assert (score >= 0) ;
	    assert (score <= n_col) ;
	    assert (head [score] >= EMPTY) ;

	    /* now add this column to dList at proper score location */
	    next_col = head [score] ;
	    Col [c].shared3.prev = EMPTY ;
	    Col [c].shared4.degree_next = next_col ;

	    /* if there already was a column with the same score, set its */
	    /* previous pointer to this new column */
	    if (next_col != EMPTY)
	    {
		Col [next_col].shared3.prev = c ;
	    }
	    head [score] = c ;

	    /* see if this score is less than current min */
	    min_score = MIN (min_score, score) ;

#ifndef NDEBUG
	    debug_count++ ;
#endif
	}
    }

#ifndef NDEBUG
    DEBUG0 (("Live cols %d out of %d, non-princ: %d\n",
	debug_count, n_col, n_col-debug_count)) ;
    assert (debug_count == n_col2) ;
    debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
#endif

    /* === Return number of remaining columns, and max row degree =========== */

    *p_n_col2 = n_col2 ;
    *p_n_row2 = n_row2 ;
    *p_max_deg = max_deg ;
}


/* ========================================================================== */
/* === find_ordering ======================================================== */
/* ========================================================================== */

/*
    Order the principal columns of the supercolumn form of the matrix
    (no supercolumns on input).  Uses a minimum approximate column minimum
    degree ordering method.  Not user-callable.
*/

PRIVATE int find_ordering	/* return the number of garbage collections */
(
    /* === Parameters ======================================================= */

    int n_row,			/* number of rows of A */
    int n_col,			/* number of columns of A */
    int Alen,			/* size of A, 2*nnz + elbow_room or larger */
    RowInfo Row [],		/* of size n_row+1 */
    ColInfo Col [],		/* of size n_col+1 */
    int A [],			/* column form and row form of A */
    int head [],		/* of size n_col+1 */
    int n_col2,			/* Remaining columns to order */
    int max_deg,		/* Maximum row degree */
    int pfree			/* index of first free slot (2*nnz on entry) */
)
{
    /* === Local variables ================================================== */

    int k ;			/* current pivot ordering step */
    int pivot_col ;		/* current pivot column */
    int *cp ;			/* a column pointer */
    int *rp ;			/* a row pointer */
    int pivot_row ;		/* current pivot row */
    int *new_cp ;		/* modified column pointer */
    int *new_rp ;		/* modified row pointer */
    int pivot_row_start ;	/* pointer to start of pivot row */
    int pivot_row_degree ;	/* # of columns in pivot row */
    int pivot_row_length ;	/* # of supercolumns in pivot row */
    int pivot_col_score ;	/* score of pivot column */
    int needed_memory ;		/* free space needed for pivot row */
    int *cp_end ;		/* pointer to the end of a column */
    int *rp_end ;		/* pointer to the end of a row */
    int row ;			/* a row index */
    int col ;			/* a column index */
    int max_score ;		/* maximum possible score */
    int cur_score ;		/* score of current column */
    unsigned int hash ;		/* hash value for supernode detection */
    int head_column ;		/* head of hash bucket */
    int first_col ;		/* first column in hash bucket */
    int tag_mark ;		/* marker value for mark array */
    int row_mark ;		/* Row [row].shared2.mark */
    int set_difference ;	/* set difference size of row with pivot row */
    int min_score ;		/* smallest column score */
    int col_thickness ;		/* "thickness" (# of columns in a supercol) */
    int max_mark ;		/* maximum value of tag_mark */
    int pivot_col_thickness ;	/* number of columns represented by pivot col */
    int prev_col ;		/* Used by Dlist operations. */
    int next_col ;		/* Used by Dlist operations. */
    int ngarbage ;		/* number of garbage collections performed */
#ifndef NDEBUG
    int debug_d ;		/* debug loop counter */
    int debug_step = 0 ;	/* debug loop counter */
#endif

    /* === Initialization and clear mark ==================================== */

    max_mark = INT_MAX - n_col ;	/* INT_MAX defined in <limits.h> */
    tag_mark = clear_mark (n_row, Row) ;
    min_score = 0 ;
    ngarbage = 0 ;
    DEBUG0 (("Ordering.. n_col2=%d\n", n_col2)) ;

    /* === Order the columns ================================================ */

    for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
    {

#ifndef NDEBUG
	if (debug_step % 100 == 0)
	{
	    DEBUG0 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
	}
	else
	{
	    DEBUG1 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
	}
	debug_step++ ;
	debug_deg_lists (n_row, n_col, Row, Col, head,
		min_score, n_col2-k, max_deg) ;
	debug_matrix (n_row, n_col, Row, Col, A) ;
#endif

	/* === Select pivot column, and order it ============================ */

	/* make sure degree list isn't empty */
	assert (min_score >= 0) ;
	assert (min_score <= n_col) ;
	assert (head [min_score] >= EMPTY) ;

#ifndef NDEBUG
	for (debug_d = 0 ; debug_d < min_score ; debug_d++)
	{
	    assert (head [debug_d] == EMPTY) ;
	}
#endif

	/* get pivot column from head of minimum degree list */
	while (head [min_score] == EMPTY && min_score < n_col)
	{
	    min_score++ ;
	}
	pivot_col = head [min_score] ;
	assert (pivot_col >= 0 && pivot_col <= n_col) ;
	next_col = Col [pivot_col].shared4.degree_next ;
	head [min_score] = next_col ;
	if (next_col != EMPTY)
	{
	    Col [next_col].shared3.prev = EMPTY ;
	}

	assert (COL_IS_ALIVE (pivot_col)) ;
	DEBUG3 (("Pivot col: %d\n", pivot_col)) ;

	/* remember score for defrag check */
	pivot_col_score = Col [pivot_col].shared2.score ;

	/* the pivot column is the kth column in the pivot order */
	Col [pivot_col].shared2.order = k ;

	/* increment order count by column thickness */
	pivot_col_thickness = Col [pivot_col].shared1.thickness ;
	k += pivot_col_thickness ;
	assert (pivot_col_thickness > 0) ;

	/* === Garbage_collection, if necessary ============================= */

	needed_memory = MIN (pivot_col_score, n_col - k) ;
	if (pfree + needed_memory >= Alen)
	{
	    pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
	    ngarbage++ ;
	    /* after garbage collection we will have enough */
	    assert (pfree + needed_memory < Alen) ;
	    /* garbage collection has wiped out the Row[].shared2.mark array */
	    tag_mark = clear_mark (n_row, Row) ;
#ifndef NDEBUG
	    debug_matrix (n_row, n_col, Row, Col, A) ;
#endif
	}

	/* === Compute pivot row pattern ==================================== */

	/* get starting location for this new merged row */
	pivot_row_start = pfree ;

	/* initialize new row counts to zero */
	pivot_row_degree = 0 ;

	/* tag pivot column as having been visited so it isn't included */
	/* in merged pivot row */
	Col [pivot_col].shared1.thickness = -pivot_col_thickness ;

	/* pivot row is the union of all rows in the pivot column pattern */
	cp = &A [Col [pivot_col].start] ;
	cp_end = cp + Col [pivot_col].length ;
	while (cp < cp_end)
	{
	    /* get a row */
	    row = *cp++ ;
	    DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
	    /* skip if row is dead */
	    if (ROW_IS_DEAD (row))
	    {
		continue ;
	    }
	    rp = &A [Row [row].start] ;
	    rp_end = rp + Row [row].length ;
	    while (rp < rp_end)
	    {
		/* get a column */
		col = *rp++ ;
		/* add the column, if alive and untagged */
		col_thickness = Col [col].shared1.thickness ;
		if (col_thickness > 0 && COL_IS_ALIVE (col))
		{
		    /* tag column in pivot row */
		    Col [col].shared1.thickness = -col_thickness ;
		    assert (pfree < Alen) ;
		    /* place column in pivot row */
		    A [pfree++] = col ;
		    pivot_row_degree += col_thickness ;
		}
	    }
	}

	/* clear tag on pivot column */
	Col [pivot_col].shared1.thickness = pivot_col_thickness ;
	max_deg = MAX (max_deg, pivot_row_degree) ;

#ifndef NDEBUG
	DEBUG3 (("check2\n")) ;
	debug_mark (n_row, Row, tag_mark, max_mark) ;
#endif

	/* === Kill all rows used to construct pivot row ==================== */

	/* also kill pivot row, temporarily */
	cp = &A [Col [pivot_col].start] ;
	cp_end = cp + Col [pivot_col].length ;
	while (cp < cp_end)
	{
	    /* may be killing an already dead row */
	    row = *cp++ ;
	    DEBUG2 (("Kill row in pivot col: %d\n", row)) ;
	    KILL_ROW (row) ;
	}

	/* === Select a row index to use as the new pivot row =============== */

	pivot_row_length = pfree - pivot_row_start ;
	if (pivot_row_length > 0)
	{
	    /* pick the "pivot" row arbitrarily (first row in col) */
	    pivot_row = A [Col [pivot_col].start] ;
	    DEBUG2 (("Pivotal row is %d\n", pivot_row)) ;
	}
	else
	{
	    /* there is no pivot row, since it is of zero length */
	    pivot_row = EMPTY ;
	    assert (pivot_row_length == 0) ;
	}
	assert (Col [pivot_col].length > 0 || pivot_row_length == 0) ;

	/* === Approximate degree computation =============================== */

	/* Here begins the computation of the approximate degree.  The column */
	/* score is the sum of the pivot row "length", plus the size of the */
	/* set differences of each row in the column minus the pattern of the */
	/* pivot row itself.  The column ("thickness") itself is also */
	/* excluded from the column score (we thus use an approximate */
	/* external degree). */

	/* The time taken by the following code (compute set differences, and */
	/* add them up) is proportional to the size of the data structure */
	/* being scanned - that is, the sum of the sizes of each column in */
	/* the pivot row.  Thus, the amortized time to compute a column score */
	/* is proportional to the size of that column (where size, in this */
	/* context, is the column "length", or the number of row indices */
	/* in that column).  The number of row indices in a column is */
	/* monotonically non-decreasing, from the length of the original */
	/* column on input to colamd. */

	/* === Compute set differences ====================================== */

	DEBUG1 (("** Computing set differences phase. **\n")) ;

	/* pivot row is currently dead - it will be revived later. */

	DEBUG2 (("Pivot row: ")) ;
	/* for each column in pivot row */
	rp = &A [pivot_row_start] ;
	rp_end = rp + pivot_row_length ;
	while (rp < rp_end)
	{
	    col = *rp++ ;
	    assert (COL_IS_ALIVE (col) && col != pivot_col) ;
	    DEBUG2 (("Col: %d\n", col)) ;

	    /* clear tags used to construct pivot row pattern */
	    col_thickness = -Col [col].shared1.thickness ;
	    assert (col_thickness > 0) ;
	    Col [col].shared1.thickness = col_thickness ;

	    /* === Remove column from degree list =========================== */

	    cur_score = Col [col].shared2.score ;
	    prev_col = Col [col].shared3.prev ;
	    next_col = Col [col].shared4.degree_next ;
	    assert (cur_score >= 0) ;
	    assert (cur_score <= n_col) ;
	    assert (cur_score >= EMPTY) ;
	    if (prev_col == EMPTY)
	    {
		head [cur_score] = next_col ;
	    }
	    else
	    {
		Col [prev_col].shared4.degree_next = next_col ;
	    }
	    if (next_col != EMPTY)
	    {
		Col [next_col].shared3.prev = prev_col ;
	    }

	    /* === Scan the column ========================================== */

	    cp = &A [Col [col].start] ;
	    cp_end = cp + Col [col].length ;