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;}statictest(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;}