old_colamd.c

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

    {
	A [i] = 0 ;
    }
    A [COLAMD_DENSE_ROW] = n_row - n_row2 ;
    A [COLAMD_DENSE_COL] = n_col - n_col2 ;
    A [COLAMD_DEFRAG_COUNT] = ngarbage ;
    A [COLAMD_JUMBLED_COLS] = init_result ;

    return (TRUE) ;
}


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

/* There are no user-callable routines beyond this point in the file */


/* ========================================================================== */
/* === init_rows_cols ======================================================= */
/* ========================================================================== */

/*
    Takes the column form of the matrix in A and creates the row form of the
    matrix.  Also, row and column attributes are stored in the Col and Row
    structs.  If the columns are un-sorted or contain duplicate row indices,
    this routine will also sort and remove duplicate row indices from the
    column form of the matrix.  Returns -1 on error, 1 if columns jumbled,
    or 0 if columns not jumbled.  Not user-callable.
*/

PRIVATE int init_rows_cols	/* returns status code */
(
    /* === Parameters ======================================================= */

    int n_row,			/* number of rows of A */
    int n_col,			/* number of columns of A */
    RowInfo Row [],		/* of size n_row+1 */
    ColInfo Col [],		/* of size n_col+1 */
    int A [],			/* row indices of A, of size Alen */
    int p []			/* pointers to columns in A, of size n_col+1 */
)
{
    /* === Local variables ================================================== */

    int col ;			/* a column index */
    int row ;			/* a row index */
    int *cp ;			/* a column pointer */
    int *cp_end ;		/* a pointer to the end of a column */
    int *rp ;			/* a row pointer */
    int *rp_end ;		/* a pointer to the end of a row */
    int last_start ;		/* start index of previous column in A */
    int start ;			/* start index of column in A */
    int last_row ;		/* previous row */
    int jumbled_columns ;	/* indicates if columns are jumbled */

    /* === Initialize columns, and check column pointers ==================== */

    last_start = 0 ;
    for (col = 0 ; col < n_col ; col++)
    {
	start = p [col] ;
	if (start < last_start)
	{
	    /* column pointers must be non-decreasing */
	    DEBUG0 (("colamd error!  last p %d p [col] %d\n",last_start,start));
	    return (-1) ;
	}
	Col [col].start = start ;
	Col [col].length = p [col+1] - start ;
	Col [col].shared1.thickness = 1 ;
	Col [col].shared2.score = 0 ;
	Col [col].shared3.prev = EMPTY ;
	Col [col].shared4.degree_next = EMPTY ;
	last_start = start ;
    }
    /* must check the end pointer for last column */
    if (p [n_col] < last_start)
    {
	/* column pointers must be non-decreasing */
	DEBUG0 (("colamd error!  last p %d p [n_col] %d\n",p[col],last_start)) ;
	return (-1) ;
    }

    /* p [0..n_col] no longer needed, used as "head" in subsequent routines */

    /* === Scan columns, compute row degrees, and check row indices ========= */

    jumbled_columns = FALSE ;

    for (row = 0 ; row < n_row ; row++)
    {
	Row [row].length = 0 ;
	Row [row].shared2.mark = -1 ;
    }

    for (col = 0 ; col < n_col ; col++)
    {
	last_row = -1 ;

	cp = &A [p [col]] ;
	cp_end = &A [p [col+1]] ;

	while (cp < cp_end)
	{
	    row = *cp++ ;

	    /* make sure row indices within range */
	    if (row < 0 || row >= n_row)
	    {
		DEBUG0 (("colamd error!  col %d row %d last_row %d\n",
			 col, row, last_row)) ;
		return (-1) ;
	    }
	    else if (row <= last_row)
	    {
		/* row indices are not sorted or repeated, thus cols */
		/* are jumbled */
		jumbled_columns = TRUE ;
	    }
	    /* prevent repeated row from being counted */
	    if (Row [row].shared2.mark != col)
	    {
		Row [row].length++ ;
		Row [row].shared2.mark = col ;
		last_row = row ;
	    }
	    else
	    {
		/* this is a repeated entry in the column, */
		/* it will be removed */
		Col [col].length-- ;
	    }
	}
    }

    /* === Compute row pointers ============================================= */

    /* row form of the matrix starts directly after the column */
    /* form of matrix in A */
    Row [0].start = p [n_col] ;
    Row [0].shared1.p = Row [0].start ;
    Row [0].shared2.mark = -1 ;
    for (row = 1 ; row < n_row ; row++)
    {
	Row [row].start = Row [row-1].start + Row [row-1].length ;
	Row [row].shared1.p = Row [row].start ;
	Row [row].shared2.mark = -1 ;
    }

    /* === Create row form ================================================== */

    if (jumbled_columns)
    {
	/* if cols jumbled, watch for repeated row indices */
	for (col = 0 ; col < n_col ; col++)
	{
	    cp = &A [p [col]] ;
	    cp_end = &A [p [col+1]] ;
	    while (cp < cp_end)
	    {
		row = *cp++ ;
		if (Row [row].shared2.mark != col)
		{
		    A [(Row [row].shared1.p)++] = col ;
		    Row [row].shared2.mark = col ;
		}
	    }
	}
    }
    else
    {
	/* if cols not jumbled, we don't need the mark (this is faster) */
	for (col = 0 ; col < n_col ; col++)
	{
	    cp = &A [p [col]] ;
	    cp_end = &A [p [col+1]] ;
	    while (cp < cp_end)
	    {
		A [(Row [*cp++].shared1.p)++] = col ;
	    }
	}
    }

    /* === Clear the row marks and set row degrees ========================== */

    for (row = 0 ; row < n_row ; row++)
    {
	Row [row].shared2.mark = 0 ;
	Row [row].shared1.degree = Row [row].length ;
    }

    /* === See if we need to re-create columns ============================== */

    if (jumbled_columns)
    {

#ifndef NDEBUG
	/* make sure column lengths are correct */
	for (col = 0 ; col < n_col ; col++)
	{
	    p [col] = Col [col].length ;
	}
	for (row = 0 ; row < n_row ; row++)
	{
	    rp = &A [Row [row].start] ;
	    rp_end = rp + Row [row].length ;
	    while (rp < rp_end)
	    {
		p [*rp++]-- ;
	    }
	}
	for (col = 0 ; col < n_col ; col++)
	{
	    assert (p [col] == 0) ;
	}
	/* now p is all zero (different than when debugging is turned off) */
#endif

	/* === Compute col pointers ========================================= */

	/* col form of the matrix starts at A [0]. */
	/* Note, we may have a gap between the col form and the row */
	/* form if there were duplicate entries, if so, it will be */
	/* removed upon the first garbage collection */
	Col [0].start = 0 ;
	p [0] = Col [0].start ;
	for (col = 1 ; col < n_col ; col++)
	{
	    /* note that the lengths here are for pruned columns, i.e. */
	    /* no duplicate row indices will exist for these columns */
	    Col [col].start = Col [col-1].start + Col [col-1].length ;
	    p [col] = Col [col].start ;
	}

	/* === Re-create col form =========================================== */

	for (row = 0 ; row < n_row ; row++)
	{
	    rp = &A [Row [row].start] ;
	    rp_end = rp + Row [row].length ;
	    while (rp < rp_end)
	    {
		A [(p [*rp++])++] = row ;
	    }
	}
	return (1) ;
    }
    else
    {
	/* no columns jumbled (this is faster) */
	return (0) ;
    }
}


/* ========================================================================== */
/* === init_scoring ========================================================= */
/* ========================================================================== */

/*
    Kills dense or empty columns and rows, calculates an initial score for
    each column, and places all columns in the degree lists.  Not user-callable.
*/

PRIVATE void init_scoring
(
    /* === Parameters ======================================================= */

    int n_row,			/* number of rows of A */
    int n_col,			/* number of columns of A */
    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 */
    double knobs [COLAMD_KNOBS],/* parameters */
    int *p_n_row2,		/* number of non-dense, non-empty rows */
    int *p_n_col2,		/* number of non-dense, non-empty columns */
    int *p_max_deg		/* maximum row degree */
)
{
    /* === Local variables ================================================== */

    int c ;			/* a column index */
    int r, row ;		/* a row index */
    int *cp ;			/* a column pointer */
    int deg ;			/* degree (# entries) of a row or column */
    int *cp_end ;		/* a pointer to the end of a column */
    int *new_cp ;		/* new column pointer */
    int col_length ;		/* length of pruned column */
    int score ;			/* current column score */
    int n_col2 ;		/* number of non-dense, non-empty columns */
    int n_row2 ;		/* number of non-dense, non-empty rows */
    int dense_row_count ;	/* remove rows with more entries than this */
    int dense_col_count ;	/* remove cols with more entries than this */
    int min_score ;		/* smallest column score */
    int max_deg ;		/* maximum row degree */
    int next_col ;		/* Used to add to degree list.*/
#ifndef NDEBUG
    int debug_count ;		/* debug only. */
#endif

    /* === Extract knobs ==================================================== */

    dense_row_count = MAX (0, MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ;
    dense_col_count = MAX (0, MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ;
    DEBUG0 (("densecount: %d %d\n", dense_row_count, dense_col_count)) ;
    max_deg = 0 ;
    n_col2 = n_col ;
    n_row2 = n_row ;

    /* === Kill empty columns =============================================== */

    /* Put the empty columns at the end in their natural, so that LU */
    /* factorization can proceed as far as possible. */
    for (c = n_col-1 ; c >= 0 ; c--)
    {
	deg = Col [c].length ;
	if (deg == 0)
	{
	    /* this is a empty column, kill and order it last */
	    Col [c].shared2.order = --n_col2 ;
	    KILL_PRINCIPAL_COL (c) ;
	}
    }
    DEBUG0 (("null columns killed: %d\n", n_col - n_col2)) ;

    /* === Kill dense columns =============================================== */

    /* Put the dense columns at the end, in their natural order */
    for (c = n_col-1 ; c >= 0 ; c--)
    {
	/* skip any dead columns */
	if (COL_IS_DEAD (c))
	{
	    continue ;
	}
	deg = Col [c].length ;
	if (deg > dense_col_count)
	{
	    /* this is a dense column, kill and order it last */
	    Col [c].shared2.order = --n_col2 ;
	    /* decrement the row degrees */
	    cp = &A [Col [c].start] ;
	    cp_end = cp + Col [c].length ;
	    while (cp < cp_end)
	    {
		Row [*cp++].shared1.degree-- ;
	    }
	    KILL_PRINCIPAL_COL (c) ;
	}
    }
    DEBUG0 (("Dense and null columns killed: %d\n", n_col - n_col2)) ;

    /* === Kill dense and empty rows ======================================== */

    for (r = 0 ; r < n_row ; r++)
    {
	deg = Row [r].shared1.degree ;
	assert (deg >= 0 && deg <= n_col) ;
	if (deg > dense_row_count || deg == 0)
	{
	    /* kill a dense or empty row */
	    KILL_ROW (r) ;
	    --n_row2 ;
	}
	else
	{
	    /* keep track of max degree of remaining rows */
	    max_deg = MAX (max_deg, deg) ;
	}
    }
    DEBUG0 (("Dense and null rows killed: %d\n", n_row - n_row2)) ;

    /* === Compute initial column scores ==================================== */

    /* At this point the row degrees are accurate.  They reflect the number */
    /* of "live" (non-dense) columns in each row.  No empty rows exist. */
    /* Some "live" columns may contain only dead rows, however.  These are */
    /* pruned in the code below. */

    /* now find the initial matlab score for each column */
    for (c = n_col-1 ; c >= 0 ; c--)
    {
	/* skip dead column */
	if (COL_IS_DEAD (c))
	{
	    continue ;
	}
	score = 0 ;
	cp = &A [Col [c].start] ;
	new_cp = cp ;
	cp_end = cp + Col [c].length ;
	while (cp < cp_end)
	{
	    /* get a row */
	    row = *cp++ ;
	    /* skip if dead */
	    if (ROW_IS_DEAD (row))
	    {
		continue ;
	    }
	    /* compact the column */