spbkp.c

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

    if ( i < 0 || i > A->n || j < 0 || j >= A->n )
	error(E_BOUNDS,"max_row_col");

    max_val = 0.0;

    idx = unord_get_idx(&(A->row[i]),j);
    if ( idx < 0 )
    {
	row_num = -1;	idx = j;
	e = chase_past(A,j,&row_num,&idx,i);
    }
    else
    {
	row_num = i;
	e = &(A->row[i].elt[idx]);
    }
    while ( row_num >= 0 && row_num < j )
    {
	if ( row_num != l )
	{
	    tmp = fabs(e->val);
	    if ( tmp > max_val )
		max_val = tmp;
	}
	e = bump_col(A,j,&row_num,&idx);
    }
    r = &(A->row[j]);
    for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
    {
	if ( e->col > j && e->col != l )
	{
	    tmp = fabs(e->val);
	    if ( tmp > max_val )
		max_val = tmp;
	}
    }

    return max_val;
}

/* nonzeros -- counts non-zeros in A */
#ifndef ANSI_C
static int	nonzeros(A)
SPMAT	*A;
#else
static int	nonzeros(const SPMAT *A)
#endif
{
    int		cnt, i;

    if ( ! A )
	return 0;
    cnt = 0;
    for ( i = 0; i < A->m; i++ )
	cnt += A->row[i].len;

    return cnt;
}

/* chk_col_access -- for spBKPfactor()
	-- checks that column access path is OK */
#ifndef ANSI_C
int	chk_col_access(A)
SPMAT	*A;
#else
int	chk_col_access(const SPMAT *A)
#endif
{
    int		cnt_nz, j, row, idx;
    SPROW	*r;
    row_elt	*e;

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

    /* count nonzeros as we go down columns */
    cnt_nz = 0;
    for ( j = 0; j < A->n; j++ )
    {
	row = A->start_row[j];
	idx = A->start_idx[j];
	while ( row >= 0 )
	{
	    if ( row >= A->m || idx < 0 )
		return FALSE;
	    r = &(A->row[row]);
	    if ( idx >= r->len )
		return FALSE;
	    e = &(r->elt[idx]);
	    if ( e->nxt_row >= 0 && e->nxt_row <= row )
		return FALSE;
	    row = e->nxt_row;
	    idx = e->nxt_idx;
	    cnt_nz++;
	}
    }

    if ( cnt_nz != nonzeros(A) )
	return FALSE;
    else
	return TRUE;
}

/* col_cmp -- compare two columns -- for sorting rows using qsort() */
#ifndef ANSI_C
static int	col_cmp(e1,e2)
row_elt	*e1, *e2;
#else
static int	col_cmp(const row_elt *e1, const row_elt *e2)
#endif
{
    return e1->col - e2->col;
}

/* spBKPfactor -- sparse Bunch-Kaufman-Parlett factorisation of A in-situ
   -- A is factored into the form P'AP = MDM' where 
   P is a permutation matrix, M lower triangular and D is block
   diagonal with blocks of size 1 or 2
   -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */
#ifndef ANSI_C
SPMAT	*spBKPfactor(A,pivot,blocks,tol)
SPMAT	*A;
PERM	*pivot, *blocks;
double	tol;
#else
SPMAT	*spBKPfactor(SPMAT *A, PERM *pivot, PERM *blocks, double tol)
#endif
{
    int		i, j, k, l, n, onebyone, r;
    int		idx, idx1, idx_piv;
    int		row_num;
    int		best_deg, best_j, best_l, best_cost, mark_cost, deg, deg_j,
			deg_l, ignore_deg;
    int		list_idx, list_idx2, old_list_idx;
    SPROW	*row, *r_piv, *r1_piv;
    row_elt	*e, *e1;
    Real	aii, aip1, aip1i;
    Real	det, max_j, max_l, s, t;
    STATIC IVEC	*scan_row = IVNULL, *scan_idx = IVNULL, *col_list = IVNULL,
		*tmp_iv = IVNULL;
    STATIC IVEC *deg_list = IVNULL;
    STATIC IVEC	*orig_idx = IVNULL, *orig1_idx = IVNULL;
    STATIC PERM	*order = PNULL;

    if ( ! A || ! pivot || ! blocks )
	error(E_NULL,"spBKPfactor");
    if ( A->m != A->n )
	error(E_SQUARE,"spBKPfactor");
    if ( A->m != pivot->size || pivot->size != blocks->size )
	error(E_SIZES,"spBKPfactor");
    if ( tol <= 0.0 || tol > 1.0 )
	error(E_RANGE,"spBKPfactor");
    
    n = A->n;
    
    px_ident(pivot);	px_ident(blocks);
    sp_col_access(A);	sp_diag_access(A);
    ignore_deg = FALSE;

    deg_list = iv_resize(deg_list,n);
    if ( order != NULL )
      px_ident(order);
    order = px_resize(order,n);
    MEM_STAT_REG(deg_list,TYPE_IVEC);
    MEM_STAT_REG(order,TYPE_PERM);

    scan_row = iv_resize(scan_row,5);
    scan_idx = iv_resize(scan_idx,5);
    col_list = iv_resize(col_list,5);
    orig_idx = iv_resize(orig_idx,5);
    orig_idx = iv_resize(orig1_idx,5);
    orig_idx = iv_resize(tmp_iv,5);
    MEM_STAT_REG(scan_row,TYPE_IVEC);
    MEM_STAT_REG(scan_idx,TYPE_IVEC);
    MEM_STAT_REG(col_list,TYPE_IVEC);
    MEM_STAT_REG(orig_idx,TYPE_IVEC);
    MEM_STAT_REG(orig1_idx,TYPE_IVEC);
    MEM_STAT_REG(tmp_iv,TYPE_IVEC);

    for ( i = 0; i < n-1; i = onebyone ? i+1 : i+2 )
    {
	/* now we want to use a Markowitz-style selection rule for
	   determining which rows to swap and whether to use
	   1x1 or 2x2 pivoting */

	/* get list of degrees of nodes */
	deg_list = iv_resize(deg_list,n-i);
	if ( ! ignore_deg )
	    for ( j = i; j < n; j++ )
		deg_list->ive[j-i] = 0;
	else
	{
	    for ( j = i; j < n; j++ )
		deg_list->ive[j-i] = 1;
	    if ( i < n )
		deg_list->ive[0] = 0;
	}
	order = px_resize(order,n-i);
	px_ident(order);

	if ( ! ignore_deg )
	{
	    for ( j = i; j < n; j++ )
	    {
		/* idx = sprow_idx(&(A->row[j]),j+1); */
		/* idx = fixindex(idx); */
		idx = 0;
		row = &(A->row[j]);
		e = &(row->elt[idx]);
		/* deg_list->ive[j-i] += row->len - idx; */
		for ( ; idx < row->len; idx++, e++ )
		    if ( e->col >= i )
			deg_list->ive[e->col - i]++;
	    }
	    /* now deg_list[k] == degree of node k+i */
	    
	    /* now sort them into increasing order */
	    iv_sort(deg_list,order);
	    /* now deg_list[idx] == degree of node i+order[idx] */
	}

	/* now we can chase through the nodes in order of increasing
	   degree, picking out the ones that satisfy our stability
	   criterion */
	list_idx = 0;	r = -1;
	best_j = best_l = -1;
	for ( deg = 0; deg <= n; deg++ )
	{
	    Real	ajj, all, ajl;

	    if ( list_idx >= deg_list->dim )
		break;	/* That's all folks! */
	    old_list_idx = list_idx;
	    while ( list_idx < deg_list->dim &&
		    deg_list->ive[list_idx] <= deg )
	    {
		j = i+order->pe[list_idx];
		if ( j < i )
		    continue;
		/* can we use row/col j for a 1 x 1 pivot? */
		/* find max_j = max_{k>=i} {|A[k][j]|,|A[j][k]|} */
		ajj = fabs(unord_get_val(A,j,j));
		if ( ajj == 0.0 )
		{
		    list_idx++;
		    continue;	/* can't use this for 1 x 1 pivot */
		}

		max_j = max_row_col(A,i,j,-1);
		if ( ajj >= tol/* *alpha */ *max_j )
		{
		    onebyone = TRUE;
		    best_j = j;
		    best_deg = deg_list->ive[list_idx];
		    break;
		}
		list_idx++;
	    }
	    if ( best_j >= 0 )
		break;
	    best_cost = 2*n;	/* > any possible Markowitz cost (bound) */
	    best_j = best_l = -1;
	    list_idx = old_list_idx;
	    while ( list_idx < deg_list->dim &&
		    deg_list->ive[list_idx] <= deg )
	    {
		j = i+order->pe[list_idx];
		ajj = fabs(unord_get_val(A,j,j));
		for ( list_idx2 = 0; list_idx2 < list_idx; list_idx2++ )
		{
		    deg_j = deg;
		    deg_l = deg_list->ive[list_idx2];
		    l = i+order->pe[list_idx2];
		    if ( l < i )
			continue;
		    /* try using rows/cols (j,l) for a 2 x 2 pivot block */
		    all = fabs(unord_get_val(A,l,l));
		    ajl = ( j > l ) ? fabs(unord_get_val(A,l,j)) :
					   fabs(unord_get_val(A,j,l));
		    det = fabs(ajj*all - ajl*ajl);
		    if ( det == 0.0 )
			continue;
		    max_j = max_row_col(A,i,j,l);
		    max_l = max_row_col(A,i,l,j);
		    if ( tol*(all*max_j+ajl*max_l) < det &&
			 tol*(ajl*max_j+ajj*max_l) < det )
		    {
			/* acceptably stable 2 x 2 pivot */
			/* this is actually an overestimate of the
			   Markowitz cost for choosing (j,l) */
			mark_cost = (ajj == 0.0) ?
			    ((all == 0.0) ? deg_j+deg_l : deg_j+2*deg_l) :
				((all == 0.0) ? 2*deg_j+deg_l :
				 2*(deg_j+deg_l));
			if ( mark_cost < best_cost )
			{
			    onebyone = FALSE;
			    best_cost = mark_cost;
			    best_j = j;
			    best_l = l;
			    best_deg = deg_j;
			}
		    }
		}
		list_idx++;
	    }
	    if ( best_j >= 0 )
		break;
	}

	if ( best_deg > (int)floor(0.8*(n-i)) )
	    ignore_deg = TRUE;

	/* now do actual interchanges */
	if ( best_j >= 0 && onebyone )
	{
	    bkp_interchange(A,i,best_j);
	    px_transp(pivot,i,best_j);
	}
	else if ( best_j >= 0 && best_l >= 0 && ! onebyone )
	{
	    if ( best_j == i || best_j == i+1 )
	    {
		if ( best_l == i || best_l == i+1 )
		{
		    /* no pivoting, but must update blocks permutation */
		    px_transp(blocks,i,i+1);
		    goto dopivot;
		}
		bkp_interchange(A,(best_j == i) ? i+1 : i,best_l);
		px_transp(pivot,(best_j == i) ? i+1 : i,best_l);
	    }
	    else if ( best_l == i || best_l == i+1 )
	    {
		bkp_interchange(A,(best_l == i) ? i+1 : i,best_j);
		px_transp(pivot,(best_l == i) ? i+1 : i,best_j);
	    }
	    else /* best_j & best_l outside i, i+1 */
	    {
		if ( i != best_j )
		{
		    bkp_interchange(A,i,best_j);
		    px_transp(pivot,i,best_j);
		}
		if ( i+1 != best_l )
		{
		    bkp_interchange(A,i+1,best_l);
		    px_transp(pivot,i+1,best_l);
		}
	    }
	}
	else	/* can't pivot &/or nothing to pivot */
	    continue;

	/* update blocks permutation */
	if ( ! onebyone )
	    px_transp(blocks,i,i+1);

	dopivot:
	if ( onebyone )
	{
	    int		idx_j, idx_k, s_idx, s_idx2;
	    row_elt	*e_ij, *e_ik;

	    r_piv = &(A->row[i]);
	    idx_piv = unord_get_idx(r_piv,i);
	    /* if idx_piv < 0 then aii == 0 and no pivoting can be done;
	       -- this means that we should continue to the next iteration */
	    if ( idx_piv < 0 )
		continue;
	    aii = r_piv->elt[idx_piv].val;
	    if ( aii == 0.0 )
		continue;

	    /* for ( j = i+1; j < n; j++ )  { ... pivot step ... } */
	    /* initialise scan_... etc for the 1 x 1 pivot */
	    scan_row = iv_resize(scan_row,r_piv->len);
	    scan_idx = iv_resize(scan_idx,r_piv->len);
	    col_list = iv_resize(col_list,r_piv->len);
	    orig_idx = iv_resize(orig_idx,r_piv->len);
	    row_num = i;	s_idx = idx = 0;
	    e = &(r_piv->elt[idx]);
	    for ( idx = 0; idx < r_piv->len; idx++, e++ )
	    {
		if ( e->col < i )
		    continue;
		scan_row->ive[s_idx] = i;
		scan_idx->ive[s_idx] = idx;
		orig_idx->ive[s_idx] = idx;
		col_list->ive[s_idx] = e->col;
		s_idx++;
	    }
	    scan_row = iv_resize(scan_row,s_idx);
	    scan_idx = iv_resize(scan_idx,s_idx);
	    col_list = iv_resize(col_list,s_idx);
	    orig_idx = iv_resize(orig_idx,s_idx);

	    order = px_resize(order,scan_row->dim);
	    px_ident(order);
	    iv_sort(col_list,order);

	    tmp_iv = iv_resize(tmp_iv,scan_row->dim);
	    for ( idx = 0; idx < order->size; idx++ )
		tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
	    iv_copy(tmp_iv,scan_idx);
	    for ( idx = 0; idx < order->size; idx++ )
		tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
	    iv_copy(tmp_iv,scan_row);
	    for ( idx = 0; idx < scan_row->dim; idx++ )
		tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
	    iv_copy(tmp_iv,orig_idx);

	    /* now do actual pivot */
	    /* for ( j = i+1; j < n-1; j++ ) .... */

	    for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
	    {
		idx_j = orig_idx->ive[s_idx];
		if ( idx_j < 0 )
		    error(E_INTERN,"spBKPfactor");
		e_ij = &(r_piv->elt[idx_j]);
		j = e_ij->col;
		if ( j < i+1 )
		    continue;
		scan_to(A,scan_row,scan_idx,col_list,j);

		/* compute multiplier */
		t = e_ij->val / aii;

		/* for ( k = j; k < n; k++ ) { .... update A[j][k] .... } */
		/* this is the row in which pivoting is done */
		row = &(A->row[j]);
		for ( s_idx2 = s_idx; s_idx2 < scan_row->dim; s_idx2++ )
		{
		    idx_k = orig_idx->ive[s_idx2];
		    e_ik = &(r_piv->elt[idx_k]);
		    k = e_ik->col;
		    /* k >= j since col_list has been sorted */

		    if ( scan_row->ive[s_idx2] == j )
		    {	/* no fill-in -- can be done directly */
			idx = scan_idx->ive[s_idx2];
			/* idx = sprow_idx2(row,k,idx); */
			row->elt[idx].val -= t*e_ik->val;
		    }
		    else
		    {	/* fill-in -- insert entry & patch column */
			int	old_row, old_idx;
			row_elt	*old_e, *new_e;

			old_row = scan_row->ive[s_idx2];
			old_idx = scan_idx->ive[s_idx2];
			/* old_idx = sprow_idx2(&(A->row[old_row]),k,old_idx); */

			if ( old_idx < 0 )
			    error(E_INTERN,"spBKPfactor");
			/* idx = sprow_idx(row,k); */
			/* idx = fixindex(idx); */
			idx = row->len;

			/* sprow_set_val(row,k,-t*e_ik->val); */
			if ( row->len >= row->maxlen )
			{ tracecatch(sprow_xpd(row,2*row->maxlen+1,TYPE_SPMAT),
				     "spBKPfactor");		}

			row->len = idx+1;

			new_e = &(row->elt[idx]);
			new_e->val = -t*e_ik->val;
			new_e->col = k;

			old_e = &(A->row[old_row].elt[old_idx]);
			new_e->nxt_row = old_e->nxt_row;
			new_e->nxt_idx = old_e->nxt_idx;
			old_e->nxt_row = j;
			old_e->nxt_idx = idx;
		    }
		}
		e_ij->val = t;
	    }
	}
	else /* onebyone == FALSE */