lhdmsym.c

speech signal process tools

lhdmsym(M1,iM1,M2,iM2,pdet,n)

/* find inverse of hermitian matrix M1+iM1, place in M2+iM2 */
/* M2 can =M1 */
/* returns 0 for non-singular matrix, 1 for singular matrix */
/* return 2 for storage allocation problems */
/* the determinant is returned as the natural log */
/* this is NOT the determinant of M1+iM1 */
/* det is zero for singular case and npd case */

register double M1[],iM1[],M2[],iM2[],*pdet;
register int n;
{
/* use ldmsyminv */

	register int i,j,iM,jM,nx2,iMx2;
	static double *Mx2;
	int rtflag;

	nx2=2*n;
	if(allocdp(&Mx2,nx2*nx2))
		{
		printf("lhdmsyminv not possible due to storage allocation problems\n"); return(2); }

/* make a 2n by 2n symmetric real matrix from the M1 and iM1 */

	for(i=iM=iMx2=0;i<n;i++,iM+=n,iMx2+=nx2)/* move M1 */
		{
		for(j=0;j<n;j++)
			Mx2[iMx2+j]=M1[iM+j];
		}
	for(i=0,iMx2=n;i<n;i++,iMx2+=nx2)/* move iM1 trans to upper cor*/
		{
		for(j=jM=0;j<n;j++,jM+=n)
			Mx2[iMx2+j]=iM1[jM+i];
		}
	for(i=iM=0,iMx2=n*nx2;i<n;i++,iM+=n,iMx2+=nx2)/* move iM1 to lower corner */
		{
		for(j=jM=0;j<n;j++)
			Mx2[iMx2+j]=iM1[iM+j];
		}
	for(i=iM=0,iMx2=n*nx2+n;i<n;i++,iM+=n,iMx2+=nx2)/* move M1 */
		{
		for(j=0;j<n;j++)
			Mx2[iMx2+j]=M1[iM+j];
		}

	rtflag=ldmsyminv(Mx2,Mx2,pdet,nx2,1);

	for(i=iM=iMx2=0;i<n;i++,iM+=n,iMx2+=nx2)
		{
		for(j=0;j<n;j++)
			M2[iM+j]=Mx2[iMx2+j];
		}
	for(i=0,iMx2=n;i<n;i++,iMx2+=nx2)
		{
		for(j=jM=0;j<n;j++,jM+=n)
			iM2[jM+i]=Mx2[iMx2+j];
		}
	freedp(&Mx2);
	return(rtflag);
}