spbkp.c

Meschach 可以解稠密或稀疏线性方程组、计算特征值和特征向量和解最小平方问题

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
  Sparse matrix Bunch--Kaufman--Parlett factorisation and solve
  Radical revision started Thu 05th Nov 1992, 09:36:12 AM
  to use Karen George's suggestion of leaving the the row elements unordered
  Radical revision completed Mon 07th Dec 1992, 10:59:57 AM
*/

static	char	rcsid[] = "$Id: spbkp.c,v 1.6 1996/08/20 19:53:10 stewart Exp $";

#include	<stdio.h>
#include	<math.h>
#include        "sparse2.h"


#ifdef MALLOCDECL
#include <malloc.h>
#endif

#define alpha	0.6403882032022076 /* = (1+sqrt(17))/8 */


#define	btos(x)	((x) ? "TRUE" : "FALSE")

/* assume no use of sqr() uses side-effects */
#define	sqr(x)	((x)*(x))

/* unord_get_idx -- returns index (encoded if entry not allocated)
	of the element of row r with column j
	-- uses linear search */
#ifndef ANSI_C
int	unord_get_idx(r,j)
SPROW	*r;
int	j;
#else
int	unord_get_idx(SPROW *r, int j)
#endif
{
    int		idx;
    row_elt	*e;

    if ( ! r || ! r->elt )
	error(E_NULL,"unord_get_idx");
    for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
	if ( e->col == j )
	    break;
    if ( idx >= r->len )
	return -(r->len+2);
    else
	return idx;
}

/* unord_get_val -- returns value of the (i,j) entry of A
	-- same assumptions as unord_get_idx() */
#ifndef ANSI_C
double	unord_get_val(A,i,j)
SPMAT	*A;
int	i, j;
#else
double	unord_get_val(SPMAT *A, int i, int j)
#endif
{
    SPROW	*r;
    int		idx;

    if ( ! A )
	error(E_NULL,"unord_get_val");
    if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
	error(E_BOUNDS,"unord_get_val");

    r = &(A->row[i]);
    idx = unord_get_idx(r,j);
    if ( idx < 0 )
	return 0.0;
    else
	return r->elt[idx].val;
}

	    
/* bkp_swap_elt -- swaps the (i,j) with the (k,l) entry of sparse matrix
	-- either or both of the entries may be unallocated */
#ifndef ANSI_C
static SPMAT	*bkp_swap_elt(A,i1,j1,idx1,i2,j2,idx2)
SPMAT	*A;
int	i1, j1, idx1, i2, j2, idx2;
#else
static SPMAT	*bkp_swap_elt(SPMAT *A, int i1, int j1, 
			      int idx1, int i2, int j2, int idx2)
#endif
{
    int		tmp_row, tmp_idx;
    SPROW	*r1, *r2;
    row_elt	*e1, *e2;
    Real	tmp;

    if ( ! A )
	error(E_NULL,"bkp_swap_elt");

    if ( i1 < 0 || j1 < 0 || i2 < 0 || j2 < 0 ||
	 i1 >= A->m || j1 >= A->n || i2 >= A->m || j2 >= A->n )
    {
	error(E_BOUNDS,"bkp_swap_elt");
    }

    if ( i1 == i2 && j1 == j2 )
	return A;
    if ( idx1 < 0 && idx2 < 0 )	/* neither allocated */
	return A;

    r1 = &(A->row[i1]);		r2 = &(A->row[i2]);
    /* if ( idx1 >= r1->len || idx2 >= r2->len )
	error(E_BOUNDS,"bkp_swap_elt"); */
    if ( idx1 < 0 )	/* assume not allocated */
    {
	idx1 = r1->len;
	if ( idx1 >= r1->maxlen )
	{    tracecatch(sprow_xpd(r1,2*r1->maxlen+1,TYPE_SPMAT),
			"bkp_swap_elt");	}
	r1->len = idx1+1;
	r1->elt[idx1].col = j1;
	r1->elt[idx1].val = 0.0;
	/* now patch up column access path */
	tmp_row = -1;	tmp_idx = j1;
	chase_col(A,j1,&tmp_row,&tmp_idx,i1-1);

	if ( tmp_row < 0 )
	{
	    r1->elt[idx1].nxt_row = A->start_row[j1];
	    r1->elt[idx1].nxt_idx = A->start_idx[j1];
	    A->start_row[j1] = i1;
	    A->start_idx[j1] = idx1;
	}
	else
	{
	    row_elt	*tmp_e;

	    tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
	    r1->elt[idx1].nxt_row = tmp_e->nxt_row;
	    r1->elt[idx1].nxt_idx = tmp_e->nxt_idx;
	    tmp_e->nxt_row = i1;
	    tmp_e->nxt_idx = idx1;
	}
    }
    else if ( r1->elt[idx1].col != j1 )
	error(E_INTERN,"bkp_swap_elt");
    if ( idx2 < 0 )
    {
	idx2 = r2->len;
	if ( idx2 >= r2->maxlen )
	{    tracecatch(sprow_xpd(r2,2*r2->maxlen+1,TYPE_SPMAT),
			"bkp_swap_elt");	}

	r2->len = idx2+1;
	r2->elt[idx2].col = j2;
	r2->elt[idx2].val = 0.0;
	/* now patch up column access path */
	tmp_row = -1;	tmp_idx = j2;
	chase_col(A,j2,&tmp_row,&tmp_idx,i2-1);
	if ( tmp_row < 0 )
	{
	    r2->elt[idx2].nxt_row = A->start_row[j2];
	    r2->elt[idx2].nxt_idx = A->start_idx[j2];
	    A->start_row[j2] = i2;
	    A->start_idx[j2] = idx2;
	}
	else
	{
	    row_elt	*tmp_e;

	    tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
	    r2->elt[idx2].nxt_row = tmp_e->nxt_row;
	    r2->elt[idx2].nxt_idx = tmp_e->nxt_idx;
	    tmp_e->nxt_row = i2;
	    tmp_e->nxt_idx = idx2;
	}
    }
    else if ( r2->elt[idx2].col != j2 )
	error(E_INTERN,"bkp_swap_elt");

    e1 = &(r1->elt[idx1]);	e2 = &(r2->elt[idx2]);

    tmp = e1->val;
    e1->val = e2->val;
    e2->val = tmp;

    return A;
}

/* bkp_bump_col -- bumps row and idx to next entry in column j */
#ifndef ANSI_C
row_elt	*bkp_bump_col(A, j, row, idx)
SPMAT	*A;
int	j, *row, *idx;
#else
row_elt	*bkp_bump_col(SPMAT *A, int j, int *row, int *idx)
#endif
{
    SPROW	*r;
    row_elt	*e;

    if ( *row < 0 )
    {
	*row = A->start_row[j];
	*idx = A->start_idx[j];
    }
    else
    {
	r = &(A->row[*row]);
	e = &(r->elt[*idx]);
	if ( e->col != j )
	    error(E_INTERN,"bkp_bump_col");
	*row = e->nxt_row;
	*idx = e->nxt_idx;
    }
    if ( *row < 0 )
	return (row_elt *)NULL;
    else
	return &(A->row[*row].elt[*idx]);
}

/* bkp_interchange -- swap rows/cols i and j (symmetric pivot)
	-- uses just the upper triangular part */
#ifndef ANSI_C
SPMAT	*bkp_interchange(A, i1, i2)
SPMAT	*A;
int	i1, i2;
#else
SPMAT	*bkp_interchange(SPMAT *A, int i1, int i2)
#endif
{
    int		tmp_row, tmp_idx;
    int		row1, row2, idx1, idx2, tmp_row1, tmp_idx1, tmp_row2, tmp_idx2;
    SPROW	*r1, *r2;
    row_elt	*e1, *e2;
    IVEC	*done_list = IVNULL;

    if ( ! A )
	error(E_NULL,"bkp_interchange");
    if ( i1 < 0 || i1 >= A->n || i2 < 0 || i2 >= A->n )
	error(E_BOUNDS,"bkp_interchange");
    if ( A->m != A->n )
	error(E_SQUARE,"bkp_interchange");

    if ( i1 == i2 )
	return A;
    if ( i1 > i2 )
    {	tmp_idx = i1;	i1 = i2;	i2 = tmp_idx;	}

    done_list = iv_resize(done_list,A->n);
    for ( tmp_idx = 0; tmp_idx < A->n; tmp_idx++ )
	done_list->ive[tmp_idx] = FALSE;
    row1 = -1;		idx1 = i1;
    row2 = -1;		idx2 = i2;
    e1 = bkp_bump_col(A,i1,&row1,&idx1);
    e2 = bkp_bump_col(A,i2,&row2,&idx2);

    while ( (row1 >= 0 && row1 < i1) || (row2 >= 0 && row2 < i1) )
	/* Note: "row2 < i1" not "row2 < i2" as we must stop before the
	   "knee bend" */
    {
	if ( row1 >= 0 && row1 < i1 && ( row1 < row2 || row2 < 0 ) )
	{
	    tmp_row1 = row1;	tmp_idx1 = idx1;
	    e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
	    if ( ! done_list->ive[row1] )
	    {
		if ( row1 == row2 )
		    bkp_swap_elt(A,row1,i1,idx1,row1,i2,idx2);
		else
		    bkp_swap_elt(A,row1,i1,idx1,row1,i2,-1);
		done_list->ive[row1] = TRUE;
	    }
	    row1 = tmp_row1;	idx1 = tmp_idx1;
	}
	else if ( row2 >= 0 && row2 < i1 && ( row2 < row1 || row1 < 0 ) )
	{
	    tmp_row2 = row2;	tmp_idx2 = idx2;
	    e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
	    if ( ! done_list->ive[row2] )
	    {
		if ( row1 == row2 )
		    bkp_swap_elt(A,row2,i1,idx1,row2,i2,idx2);
		else
		    bkp_swap_elt(A,row2,i1,-1,row2,i2,idx2);
		done_list->ive[row2] = TRUE;
	    }
	    row2 = tmp_row2;	idx2 = tmp_idx2;
	}
	else if ( row1 == row2 )
	{
	    tmp_row1 = row1;	tmp_idx1 = idx1;
	    e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
	    tmp_row2 = row2;	tmp_idx2 = idx2;
	    e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
	    if ( ! done_list->ive[row1] )
	    {
		bkp_swap_elt(A,row1,i1,idx1,row2,i2,idx2);
		done_list->ive[row1] = TRUE;
	    }
	    row1 = tmp_row1;	idx1 = tmp_idx1;
	    row2 = tmp_row2;	idx2 = tmp_idx2;
	}
    }

    /* ensure we are **past** the first knee */
    while ( row2 >= 0 && row2 <= i1 )
	e2 = bkp_bump_col(A,i2,&row2,&idx2);

    /* at/after 1st "knee bend" */
    r1 = &(A->row[i1]);
    idx1 = 0;
    e1 = &(r1->elt[idx1]);
    while ( row2 >= 0 && row2 < i2 )
    {
	/* used for update of e2 at end of loop */
	tmp_row = row2;	tmp_idx = idx2;
	if ( ! done_list->ive[row2] )
	{
	    r2 = &(A->row[row2]);
	    bkp_bump_col(A,i2,&tmp_row,&tmp_idx);
	    done_list->ive[row2] = TRUE;
	    tmp_idx1 = unord_get_idx(r1,row2);
	    tracecatch(bkp_swap_elt(A,row2,i2,idx2,i1,row2,tmp_idx1),
		       "bkp_interchange");
	}

	/* update e1 and e2 */
	row2 = tmp_row;	idx2 = tmp_idx;
	e2 = ( row2 >= 0 ) ? &(A->row[row2].elt[idx2]) : (row_elt *)NULL;
    }

    idx1 = 0;
    e1 = r1->elt;
    while ( idx1 < r1->len )
    {
	if ( e1->col >= i2 || e1->col <= i1 )
	{
	    idx1++;
	    e1++;
	    continue;
	}
	if ( ! done_list->ive[e1->col] )
	{
	    tmp_idx2 = unord_get_idx(&(A->row[e1->col]),i2);
	    tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,e1->col,i2,tmp_idx2),
		       "bkp_interchange");
	    done_list->ive[e1->col] = TRUE;
	}
	idx1++;
	e1++;
    }

    /* at/after 2nd "knee bend" */
    idx1 = 0;
    e1 = &(r1->elt[idx1]);
    r2 = &(A->row[i2]);
    idx2 = 0;
    e2 = &(r2->elt[idx2]);
    while ( idx1 < r1->len )
    {
	if ( e1->col <= i2 )
	{
	    idx1++;	e1++;
	    continue;
	}
	if ( ! done_list->ive[e1->col] )
	{
	    tmp_idx2 = unord_get_idx(r2,e1->col);
	    tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,i2,e1->col,tmp_idx2),
		       "bkp_interchange");
	    done_list->ive[e1->col] = TRUE;
	}
	idx1++;	e1++;
    }

    idx2 = 0;	e2 = r2->elt;
    while ( idx2 < r2->len )
    {
	if ( e2->col <= i2 )
	{
	    idx2++;	e2++;
	    continue;
	}
	if ( ! done_list->ive[e2->col] )
	{
	    tmp_idx1 = unord_get_idx(r1,e2->col);
	    tracecatch(bkp_swap_elt(A,i2,e2->col,idx2,i1,e2->col,tmp_idx1),
		       "bkp_interchange");
	    done_list->ive[e2->col] = TRUE;
	}
	idx2++;	e2++;
    }

    /* now interchange the digonal entries! */
    idx1 = unord_get_idx(&(A->row[i1]),i1);
    idx2 = unord_get_idx(&(A->row[i2]),i2);
    if ( idx1 >= 0 || idx2 >= 0 )
    {
	tracecatch(bkp_swap_elt(A,i1,i1,idx1,i2,i2,idx2),
		   "bkp_interchange");
    }

    return A;
}


/* iv_min -- returns minimum of an integer vector
   -- sets index to the position in iv if index != NULL */
#ifndef ANSI_C
int	iv_min(iv,index)
IVEC	*iv;
int	*index;
#else
int	iv_min(IVEC *iv, int *index)
#endif
{
    int		i, i_min, min_val, tmp;
    
    if ( ! iv ) 
	error(E_NULL,"iv_min");
    if ( iv->dim <= 0 )
	error(E_SIZES,"iv_min");
    i_min = 0;
    min_val = iv->ive[0];
    for ( i = 1; i < iv->dim; i++ )
    {
	tmp = iv->ive[i];
	if ( tmp < min_val )
	{
	    min_val = tmp;
	    i_min = i;
	}
    }
    
    if ( index != (int *)NULL )
	*index = i_min;
    
    return min_val;
}

/* max_row_col -- returns max { |A[j][k]| : k >= i, k != j, k != l } given j
	using symmetry and only the upper triangular part of A */
#ifndef ANSI_C
static double max_row_col(A,i,j,l)
SPMAT	*A;
int	i, j, l;
#else
static double max_row_col(SPMAT *A, int i,int j, int l)
#endif
{
    int		row_num, idx;
    SPROW	*r;
    row_elt	*e;
    Real	max_val, tmp;

    if ( ! A )
	error(E_NULL,"max_row_col");