lsyminv.c

speech signal process tools

lsyminv(M,Minv,Temp,n,pdet)	/* M symmetric (toeplitz?) positive definite */
				/* Minv=Ainv*Ainvt */
		/* cholesky finds A, such that M=At*A, A upper triangular */
		/* triinv finds Ainv */
		/* return 0 normal, 1 M singular, 2 M negative definite (???) */
		/* returns ln det of M in pdet */
register int n;
double M[],Minv[],Temp[],*pdet;
{
	register double sum;
	register int i,j,k,iM,jM;
	int rtflag;
	double log();

	for(i=iM=0;i<n;i++,iM+=n)
		{
		for(j=0;j<n;j++)
			Temp[iM+j]=M[iM+j];
		}
	if((rtflag=cholesky(Temp,Minv,n)))
		{
		*pdet=0;
		return(rtflag);
		}
	/* Minv contains A, Temp scribbled */
	if(triinv(Minv,Temp,n))
		return(1);
		/* Temp==Ainv, Minv==A not scribbled, do det */
	sum=0;
	for(i=iM=0;i<n;i++,iM+=n)
		sum+=log(Minv[iM+i]);
	*pdet=2*sum;/* return ln det M */
	for(i=iM=0;i<n;i++,iM+=n)
		{
		for(j=i,jM=i*n;j<n;j++,jM+=n)
			{
			sum=0;
			for(k=0;k<n;k++)
				sum+=Temp[iM+k]*Temp[jM+k];
			Minv[jM+i]=Minv[iM+j]=sum;
			}
		}
	return(0);
}