ldmsyminv.c

speech signal process tools

ldmsyminv(M1,M2,pdet,n,flag)
/* this routine appears to give the inverse for non symmetric matrices */
/* produces the inverse of the symmetric matrix M1 in M2 */
/* M2 can =M1, M1 is destroyed unless flag=1 */
/* returns 0 for non-singular matrix, 1 for singular matrix */
/* return 2 for storage allocation problems */
/* the determinant is returned as the natural log */
/* det is zero for singular case and npd case */
register double M1[],M2[],*pdet;
register int n,flag;
{
	register int rtflag,nxn;
	static double *M3,*M4;
	double det,t0,t2,t3,dtemp;
	double log();

	if(n==1)
		{
		if(M1[0]<=0)
			*pdet=0;
		else *pdet=log(M1[0]);
		if(M1[0]==0)
			return(1);
		dtemp=1/M1[0];
		M2[0]=dtemp;
		return(0);
		}
	if(n==2)
		{
		det=M1[0]*M1[3]-M1[2]*M1[1];
		if(det<=0)
			*pdet=0;
		else *pdet=log(det);
		if(det==0)
			return(1);
		t0=M1[3]/det;
		t3=M1[0]/det;
		t2= -M1[1]/det;
		M2[3]=t3;
		M2[0]=t0;
		M2[2]=M2[1]=t2;
		return(0);
		}
	nxn=n*n;
	if( (flag==1) || M1==M2 )
		{
		if(allocdp(&M3,nxn) || allocdp(&M4,nxn))
			{
			printf("ldmsyminv not possible due to storage allocation problems\n");
			return(2);
			}
		dmmove(M1,M4,n);/* move M1 to M4 */
		rtflag=lorthogin(M4,M2,M3,pdet,n);
		freedp(&M3);
		freedp(&M4);
		if(rtflag)
			return(1);
		return(0);
		}
	if(allocdp(&M3,nxn))
		{
		printf("ldmsyminv not possible due to storage allocation problems\n");
		return(2);
		}
	rtflag=lorthogin(M1,M2,M3,pdet,n);
		freedp(&M3);
	if(rtflag)
		return(1);
	return(0);
}