/* dpofa.c -- double precision solution of symmetric matrix
 *
 *	Factors the symmetric matrix A
 *
 *
 *	From "Linpack Users' Guide", Dongarra, et al.
 *	Translated from FORTRAN by C. Bond, 1994.
 */
#include <stdio.h>
#include <math.h>
#include "func.h"

#define max(a,b)	(((a) > (b)) ? (a) : (b))
#define min(a,b)	(((a) < (b)) ? (a) : (b))
 
void dpofa(double **a,int n,int *info)
{
	double t,tt,s;
	int j,jm1,k,l;
	
	for (j = 0;j < n;j++) {
		*info = j;
		s = 0.0;
		jm1 = j - 1;
		if (jm1 < 0) goto _20;
		for (k = 0; k <= jm1; k++) {
/*			ddot(k-1,a(1,k),1,a(1,j),1);	*/
			tt = 0;
			for (l = 0;l < k; l++)
				tt += (a[l][k] * a[l][j]);
				
			t = a[k][j] - tt;
			t /= a[k][k];
			a[k][j] = t;
			s += (t * t);
		}
_20:
		s = a[j][j] - s;
		if (s <= 0.0) goto _40;
		a[j][j] = sqrt(s);
	}
	*info = 0;
_40:
}