#include <math.h>
float sdot(long	n,float	*sx,long incx,float *sy,long incy)
{
static long i,ix,iy,m,mp1;
static float sdot,stemp;
    stemp = sdot = 0.0;
    if(n <= 0) return sdot;
    if(incx == 1 && incy == 1) goto S20;
    ix = iy = 1;
    if(incx < 0) ix = (-n+1)*incx+1;
    if(incy < 0) iy = (-n+1)*incy+1;
    for(i=1; i<=n; i++)	{
	stemp += (*(sx+ix-1)**(sy+iy-1));
	ix += incx;
	iy += incy;
    }
    sdot = stemp;
    return sdot;
S20:
    m =	n % 5L;
    if(m == 0) goto S40;
    for(i=0; i<m; i++) stemp +=	(*(sx+i)**(sy+i));
    if(n < 5) goto S60;
S40:
    mp1	= m+1;
    for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
      +1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
S60:
    sdot = stemp;
    return sdot;
}
void spofa(float *a,long lda,long n,long *info)
/*
     SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
     SPOFA IS USUALLY CALLED BY	SPOCO, BUT IT CAN BE CALLED
     DIRECTLY WITH A SAVING IN TIME IF	RCOND  IS NOT NEEDED.
     (TIME FOR SPOCO) =	(1 + 18/N)*(TIME FOR SPOFA) .
     ON	ENTRY
	A	REAL(LDA, N)
		THE SYMMETRIC MATRIX TO	BE FACTORED.  ONLY THE
		DIAGONAL AND UPPER TRIANGLE ARE	USED.
	LDA	INTEGER
		THE LEADING DIMENSION OF THE ARRAY  A .
	N	INTEGER
		THE ORDER OF THE MATRIX	 A .
     ON	RETURN
	A	AN UPPER TRIANGULAR MATRIX  R  SO THAT	A = TRANS(R)*R
		WHERE  TRANS(R)	 IS THE	TRANSPOSE.
		THE STRICT LOWER TRIANGLE IS UNALTERED.
		IF  INFO .NE. 0	, THE FACTORIZATION IS NOT COMPLETE.
	INFO	INTEGER
		= 0  FOR NORMAL	RETURN.
		= K  SIGNALS AN	ERROR CONDITION.  THE LEADING MINOR
		     OF	ORDER  K  IS NOT POSITIVE DEFINITE.
     LINPACK.  THIS VERSION DATED 08/14/78 .
     CLEVE MOLER, UNIVERSITY OF	NEW MEXICO, ARGONNE NATIONAL LAB.
     SUBROUTINES AND FUNCTIONS
     BLAS SDOT
     FORTRAN SQRT
     INTERNAL VARIABLES
*/
{
extern float sdot(long n,float *sx,long	incx,float *sy,long incy);
static long j,jm1,k;
static float t,s;
/*
     BEGIN BLOCK WITH ...EXITS TO 40
*/
    for(j=1; j<=n; j++)	{
	*info =	j;
	s = 0.0;
	jm1 = j-1;
	if(jm1 < 1) goto S20;
	for(k=0; k<jm1;	k++) {
	    t =	*(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
	    t /=  *(a+k+k*lda);
	    *(a+k+(j-1)*lda) = t;
	    s += (t*t);
	}
S20:
	s = *(a+j-1+(j-1)*lda)-s;
/*
     ......EXIT
*/
	if(s <=	0.0) goto S40;
	*(a+j-1+(j-1)*lda) = sqrt(s);
    }
    *info = 0;
S40:
    return;
}
