cholesky.c

speech signal process tools

cholesky(M,A,n)	/* routine destroys M but finds A such that M=At*A */
		/* M symmetric (toeplitz?) A upper triangular */
		/* 0 normal return, 1 M singular, 2 M negative definite (???) */
register int n;
register double *M,*A;
{
	register int i,j,k,iM,jM;
	double Mii,Aii;
	register double dtemp,fact;
	double sqrt();

	for(i=iM=0;i<n;i++,iM+=n)
		{
		if((Mii=M[iM+i]) == 0 )return(1);
		if(Mii<0)return(2);
		Aii=sqrt(Mii);
		A[iM+i]=Aii;
		fact=1/Aii;
		for(j=i+1;j<n;j++)
			{
			dtemp=M[iM+j]*fact;
			A[iM+j]=dtemp;
			}
		for(j=i+1,jM=iM+n;j<n;j++,jM+=n)
			{
			fact=A[iM+j];
			for(k=j;k<n;k++)
				{
				dtemp=fact*A[iM+k];
				M[jM+k]-=dtemp;
				}
			}
		}

/* clear lower triangle */

	fact=0;
	for(i=0;i<n;i++,A+=n)
		{
		for(j=0;j<i;j++)
			*(A+j)=fact;
		}
	return(0);
}