spchfctr.c

C语言版本的矩阵库

VEC	*spCHsolve(SPMAT *L, const VEC *b, VEC *out)
#endif
{
	int	i, j_idx, n, scan_idx, scan_row;
	SPROW	*row;
	row_elt	*elt;
	Real	diag_val, sum, *out_ve;

	if ( L == SMNULL || b == VNULL )
		error(E_NULL,"spCHsolve");
	if ( L->m != L->n )
		error(E_SQUARE,"spCHsolve");
	if ( b->dim != L->m )
		error(E_SIZES,"spCHsolve");

	if ( ! L->flag_col )
		sp_col_access(L);
	if ( ! L->flag_diag )
		sp_diag_access(L);

	out = v_copy(b,out);
	out_ve = out->ve;

	/* forward substitution: solve L.x=b for x */
	n = L->n;
	for ( i = 0; i < n; i++ )
	{
		sum = out_ve[i];
		row = &(L->row[i]);
		elt = row->elt;
		for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
		{
		    if ( elt->col >= i )
			break;
		    sum -= elt->val*out_ve[elt->col];
		}
		if ( row->diag >= 0 )
		    out_ve[i] = sum/(row->elt[row->diag].val);
		else
		    error(E_SING,"spCHsolve");
	}

	/* backward substitution: solve L^T.out = x for out */
	for ( i = n-1; i >= 0; i-- )
	{
		sum = out_ve[i];
		row = &(L->row[i]);
		/* Note that row->diag >= 0 by above loop */
		elt = &(row->elt[row->diag]);
		diag_val = elt->val;

		/* scan down column */
		scan_idx = elt->nxt_idx;
		scan_row = elt->nxt_row;
		while ( scan_row >= 0 /* && scan_idx >= 0 */ )
		{
		    row = &(L->row[scan_row]);
		    elt = &(row->elt[scan_idx]);
		    sum -= elt->val*out_ve[scan_row];
		    scan_idx = elt->nxt_idx;
		    scan_row = elt->nxt_row;
		}
		out_ve[i] = sum/diag_val;
	}

	return out;
}

/* spICHfactor -- sparse Incomplete Cholesky factorisation
	-- does a Cholesky factorisation assuming NO FILL-IN
	-- as for spCHfactor(), only the lower triangular part of A is used */
#ifndef ANSI_C
SPMAT	*spICHfactor(A)
SPMAT	*A;
#else
SPMAT	*spICHfactor(SPMAT *A)
#endif
{
	int	k, m, n, nxt_row, nxt_idx, diag_idx;
	Real	pivot, tmp2;
	SPROW	*r_piv, *r_op;
	row_elt	*elt_piv, *elt_op;

	if ( A == SMNULL )
		error(E_NULL,"spICHfactor");
	if ( A->m != A->n )
		error(E_SQUARE,"spICHfactor");

	/* set up access paths if not already done so */
	if ( ! A->flag_col )
		sp_col_access(A);
	if ( ! A->flag_diag )
		sp_diag_access(A);

	m = A->m;	n = A->n;
	for ( k = 0; k < m; k++ )
	{
		r_piv = &(A->row[k]);

		diag_idx = r_piv->diag;
		if ( diag_idx < 0 )
			error(E_POSDEF,"spICHfactor");

		elt_piv = r_piv->elt;

		/* set diagonal entry of Cholesky factor */
		tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
		if ( tmp2 <= 0.0 )
			error(E_POSDEF,"spICHfactor");
		elt_piv[diag_idx].val = pivot = sqrt(tmp2);

		/* find next row where something (non-trivial) happens */
		nxt_row = elt_piv[diag_idx].nxt_row;
		nxt_idx = elt_piv[diag_idx].nxt_idx;

		/* now set the k-th column of the Cholesky factors */
		while ( nxt_row >= 0 && nxt_idx >= 0 )
		{
		    /* nxt_row and nxt_idx give next next row (& index)
			of the entry to be modified */
		    r_op = &(A->row[nxt_row]);
		    elt_op = r_op->elt;
		    elt_op[nxt_idx].val = (elt_op[nxt_idx].val -
				sprow_ip(r_piv,r_op,k))/pivot;

		    nxt_row = elt_op[nxt_idx].nxt_row;
		    nxt_idx = elt_op[nxt_idx].nxt_idx;
		}
	}

	return A;
}


/* spCHsymb -- symbolic sparse Cholesky factorisation
	-- does NOT do any floating point arithmetic; just sets up the structure
	-- only the lower triangular part of A (incl. diagonal) is used */
#ifndef ANSI_C
SPMAT	*spCHsymb(A)
SPMAT	*A;
#else
SPMAT	*spCHsymb(SPMAT *A)
#endif
{
	register 	int	i;
	int	idx, k, m, minim, n, num_scan, diag_idx, tmp1;
	SPROW	*r_piv, *r_op;
	row_elt	*elt_piv, *elt_op, *old_elt;

	if ( A == SMNULL )
		error(E_NULL,"spCHsymb");
	if ( A->m != A->n )
		error(E_SQUARE,"spCHsymb");

	/* set up access paths if not already done so */
	if ( ! A->flag_col )
		sp_col_access(A);
	if ( ! A->flag_diag )
		sp_diag_access(A);

	/* printf("spCHsymb() -- checkpoint 1\n"); */
	m = A->m;	n = A->n;
	for ( k = 0; k < m; k++ )
	{
		r_piv = &(A->row[k]);
		if ( r_piv->len > scan_len )
			set_scan(r_piv->len);
		elt_piv = r_piv->elt;
		diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
		if ( diag_idx < 0 )
			error(E_POSDEF,"spCHsymb");
		old_elt = &(elt_piv[diag_idx]);
		for ( i = 0; i < r_piv->len; i++ )
		{
			if ( elt_piv[i].col > k )
				break;
			col_list[i] = elt_piv[i].col;
			scan_row[i] = elt_piv[i].nxt_row;
			scan_idx[i] = elt_piv[i].nxt_idx;
		}
		/* printf("spCHsymb() -- checkpoint 2\n"); */
		num_scan = i;	/* number of actual entries in scan_row etc. */
		/* printf("num_scan = %d\n",num_scan); */

		/* now set the k-th column of the Cholesky factors */
		/* printf("k = %d\n",k); */
		for ( ; ; )	/* forever do... */
		{
		    /* printf("spCHsymb() -- checkpoint 3\n"); */
		    /* find next row where something (non-trivial) happens
			i.e. find min(scan_row) */
		    minim = n;
		    for ( i = 0; i < num_scan; i++ )
		    {
			tmp1 = scan_row[i];
			/* printf("%d ",tmp1); */
			minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
		    }

		    if ( minim >= n )
			break;	/* nothing more to do for this column */
		    r_op = &(A->row[minim]);
		    elt_op = r_op->elt;

		    /* set next entry in column k of Cholesky factors */
		    idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
		    if ( idx < 0 )
		    {	/* fill-in */
			sp_set_val(A,minim,k,0.0);
			/* in case a realloc() has occurred... */
			elt_op = r_op->elt;
			/* now set up column access path again */
			idx = sprow_idx2(r_op,k,-(idx+2));
			tmp1 = old_elt->nxt_row;
			old_elt->nxt_row = minim;
			r_op->elt[idx].nxt_row = tmp1;
			tmp1 = old_elt->nxt_idx;
			old_elt->nxt_idx = idx;
			r_op->elt[idx].nxt_idx = tmp1;
		    }

		    /* printf("spCHsymb() -- checkpoint 4\n"); */

		    /* remember current element in column k for column chain */
		    idx = sprow_idx2(r_op,k,idx);
		    old_elt = &(r_op->elt[idx]);

		    /* update scan_row */
		    /* printf("spCHsymb() -- checkpoint 5\n"); */
		    /* printf("minim = %d\n",minim); */
		    for ( i = 0; i < num_scan; i++ )
		    {
			if ( scan_row[i] != minim )
				continue;
			idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
			if ( idx < 0 )
			{	scan_row[i] = -1;	continue;	}
			scan_row[i] = elt_op[idx].nxt_row;
			scan_idx[i] = elt_op[idx].nxt_idx;
			/* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
			/* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
		    }
			
		}
	    /* printf("spCHsymb() -- checkpoint 6\n"); */
	}

	return A;
}

/* comp_AAT -- compute A.A^T where A is a given sparse matrix */
#ifndef ANSI_C
SPMAT	*comp_AAT(A)
SPMAT	*A;
#else
SPMAT	*comp_AAT(SPMAT *A)
#endif
{
	SPMAT	*AAT;
	SPROW	*r, *r2;
	row_elt	*elts, *elts2;
	int	i, idx, idx2, j, m, minim, n, num_scan, tmp1;
	Real	ip;

	if ( ! A )
		error(E_NULL,"comp_AAT");
	m = A->m;	n = A->n;

	/* set up column access paths */
	if ( ! A->flag_col )
		sp_col_access(A);

	AAT = sp_get(m,m,10);

	for ( i = 0; i < m; i++ )
	{
		/* initialisation */
		r = &(A->row[i]);
		elts = r->elt;

		/* set up scan lists for this row */
		if ( r->len > scan_len )
		    set_scan(r->len);
		for ( j = 0; j < r->len; j++ )
		{
		    col_list[j] = elts[j].col;
		    scan_row[j] = elts[j].nxt_row;
		    scan_idx[j] = elts[j].nxt_idx;
		}
		num_scan = r->len;

		/* scan down the rows for next non-zero not
			associated with a diagonal entry */
		for ( ; ; )
		{
		    minim = m;
		    for ( idx = 0; idx < num_scan; idx++ )
		    {
			tmp1 = scan_row[idx];
			minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
		    }
		    if ( minim >= m )
		 	break;
		    r2 = &(A->row[minim]);
		    if ( minim > i )
		    {
			ip = sprow_ip(r,r2,n);
		        sp_set_val(AAT,minim,i,ip);
		        sp_set_val(AAT,i,minim,ip);
		    }
		    /* update scan entries */
		    elts2 = r2->elt;
		    for ( idx = 0; idx < num_scan; idx++ )
		    {
			if ( scan_row[idx] != minim || scan_idx[idx] < 0 )
			    continue;
			idx2 = scan_idx[idx];
			scan_row[idx] = elts2[idx2].nxt_row;
			scan_idx[idx] = elts2[idx2].nxt_idx;
		    }
		}

		/* set the diagonal entry */
		sp_set_val(AAT,i,i,sprow_sqr(r,n));
	}

	return AAT;
}