dmimprove.c

speech signal process tools

dmimprove(M,Minv,B,T,T1,n)
/* Hotelling and Bodewig method of improving on inverse of matrix */
/* see Householder page 57 */
/* M is the matrix whose inverse is desired */
/* Minv is the approximate inverse */
/* B, T,T1 are temporary matrices */
register double M[],Minv[],B[],T[],T1[];
register int n;
{
	dmmult(Minv,M,B,n);
	onedmsub(B,T,n);/* subtract B from one and put in T */
	/* T contains Householder Bo' */

	while(test(B,n))
		{
	onedmadd(T,B,n);/* add one to T, put in B */
	/* B contains 1+Bo'power */
	dmmult(B,Minv,T1,n);
	dmmove(T1,Minv,n);
	dmmove(T,B,n);
	dmmult(B,B,T,n);
	dmmult(M,Minv,B,n);
		}
	return;
}
static
test(B,n)
register double *B;
register int n;
{
	register int j,i,k;
	register double sum;
	double sqrt();

	sum=0;
	for(i=k=0;i<n;i++)
		{
		for(j=0;j<n;j++,k++)
			{
			if(i!=j)
				sum+=B[k]*B[k];
			else
				sum+=(1-B[k])*(1-B[k]);
			}
		}
	sum=sqrt(sum);
	printf("dmimprove sum=%.16e\n",sum);
	if(sum<1e-10)
		return(0);
	return(1);
}
onedmadd(B,T,n)/* add one to B, put in T */
register double *B,*T;
register int n;
{
	register int i;

	dmmove(B,T,n);
	for(i=0;i<n;i++,T+=n+1)
		*T+=1;
	return;
}
onedmsub(B,T,n)/* subtract B from one and put in T */
register double *B,*T;
register int n;
{
	register int i,j;

	for(i=0;i<n;i++)
		{
		for(j=0;j<n;j++,T++,B++)
			{
			if(i==j)
				*T= 1- *B;
			else
				*T= - *B;
			}
		}
	return;
}