 /*
  * Khoros: $Id$
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id$";
#endif

 /*
  * $Log$
  */

/*
 *   Copyright, 1991, The Regents of the University of California.
 *   This software was produced under a U.S. Government contract
 *   (W-7405-ENG-36) by the Los Alamos National Laboratory, which is
 *   operated by the University of California for the U.S. Department
 *   of Energy.  The U.S. Government is licensed to use, reproduce,
 *   and distribute this software.  Neither the Government nor the
 *   University makes any warranty, express or implied, or assumes
 *   any liability responsibility for the use of this software.
 */

/*
 *  SM_csolve ()      Solve a system of linear equations given the Cholesky
 *                    decomposition of the symmetric, positive definite
 *                    coefficient matrix A.  The matrix is assumed to
 *                    be an SM type matrix, allocated via the SM_get_matrix
 *                    routine, and NONE OF THE LOWER TRIANGULAR ELEMENTS
 *                    are ever modified.
 *
 *  WARNING:  NO ERROR CHECKING IS PERFORMED BY THIS ROUTINE AT THE MOMENT!
 *
 *  Written by:  Patrick M. Kelly
 *  Date:        2/25/91
 */

void SM_csolve ( dim, R, b )

int dim;			/* Dimension of the problem */
float **R;		/* Upper triangular matrix of Cholesky factorization */
float *b;			/* (IN) Right hand side   (OUT) Solution vector */

{
	int i, j;		/* Loop control */

	/*
	 *      Solve Ax = b        A = transpose (R) * R
	 */

	/*
	 *  FORWARD SUBSTITUTION      Ly = b    L = transpose (R)
	 */
	for ( i = 0 ; i < dim ; i ++ ) {
		for ( j = 0 ; j < i ; j ++ ) {
			b [i] -= R [j][i] * b [j];
		}
		b [i] /= R [i][i];
	}

	/*
	 *  BACKWARD SUBSTITUTION      Ux = y   U = R
	 */
	for ( i = dim-1 ; i >= 0 ; i -- ) {
		for ( j = i+1 ; j < dim ; j ++ ) {
			b [i] -= R [i][j] * b [j];
		}
		b [i] /= R [i][i];
	}
	
	return;
}
