lorthogin.c

来自「speech signal process tools」· C语言 代码 · 共 47 行

lorthogin(M,Minv,Tinv,pdet,n)
register double *M,*Minv,*Tinv,*pdet;
register int n;
/* returns log(det) , 0 if negative definite or singular */
{
	register int i,j,toggle;
	register double Mt,*s,*s1;
	double log(),dot();

	if(orthog(M,Minv,n))
		{
		*pdet=0;
		return(1);
		}
	/* orthog returns an orthogonal matrix O in M and an upper triagonal
	matrix T in Minv such that M=O*T */

	*pdet=0;
	toggle=1;
	j=n+1;
	for(i=0,s=Minv;i<n;i++,s+=j)
		{
		if((Mt= *s)<0)
			{
			toggle= -toggle;
			Mt= -Mt;
			}
		if(!Mt)
			{
			toggle= -1;
			break;
			}
		*pdet+=log(Mt);
		}
	if(toggle!=1)
		*pdet=0;

	if(triinv(Minv,Tinv,n))
		return(1);
	for(i=0,s=Minv;i<n;i++,Tinv+=n)
		{
		for(j=0,s1=M;j<n;j++,s1+=n,s++)
			*s=dot(Tinv,1,s1,1,n);
		}
	return(0);
}