 /*
  * 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.
 */

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

/*
 *  SM_cfactor ()     Calculate the Cholesky factorization of a symmetric,
 *                    positive definite matrix.  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.
 *
 *                    Upon routine exit, the matrix A is assumed to contain
 *                    an UPPER TRIANGULAR MATRIX R (the zero elements below the
 *                    diagonal SHOULD NEVER BE INDEXED due to the structure
 *                    of SM type matrices) such that:
 *
 *                                 A = transpose (R) * R
 *
 *  WARNING:  LITTLE ERROR CHECKING IS PERFORMED BY THIS ROUTINE AT THE MOMENT!
 *            INPUT MATRIX A *MUST* BE POSITIVE DEFINITE UPON ENTRY.
 *
 *  Written by:  Patrick M. Kelly
 *  Date:        2/25/91
 *
 *  MODIFICATIONS:
 *
 *     7/31/91  Patrick M. Kelly
 *
 *              Changed routine to return an integer instead of a void
 *              routine.  Now returns a 1 upon successful factorization, 
 *              and 0 otherwise.
 */

int SM_cfactor ( dim, A )

int dim;				/* Dimension of the problem */
float **A;			/* The SM matrix to factor */

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

	for ( k = 0 ; k < dim ; k ++ ) {

		if ( A [k][k] <= 0.0 ) return ( 0 );

		A [k][k] = sqrt ( A [k][k] ) ;
		for ( i = k+1 ; i < dim ; i ++ ) {
			A [k][i] = A [k][i] / A [k][k] ;
			for ( j = k+1 ; j <= i ; j ++ ) {
				A [j][i] -= A [k][i] * A [k][j] ;
			}
		}

	}

	return ( 1 );
}
