matrixinverse.c

C语言实现复数型矩阵求逆

#include "MatrixInverse.h"

/*float complex vector */

floatComplexMatrix *floatComplexMatrixAlloc(size_t row,size_t col)
{
  size_t i;
  floatComplexMatrix *r;
  r=malloc(sizeof(floatComplexMatrix));
  r->a=malloc(row*sizeof(floatComplex *));
  for(i=0;i<row;i++)
    r->a[i]=malloc(col*sizeof(floatComplex));
  r->row=row;
  r->col=col;
  return r;
}

void floatComplexMatrixFree(floatComplexMatrix *mat)
{
  size_t i;
  for(i=0;i<mat->row;i++)
    free(mat->a[i]);
  free(mat->a);
  free(mat);
}
void floatPrintComplexMatrix(floatComplexMatrix *mat)
{
	int i,j;
	for(i=0;i<mat->row;i++)
	{
		for(j=0;j<mat->col;j++)
		{
			printf("%f+i(%f)  ",mat->a[i][j].re,mat->a[i][j].im);
		}
		printf("\n");
	}
}
/* float complex number */

floatComplex floatComplexSet(float real,float imag)
{
  floatComplex a;
  a.re=real;
  a.im=imag;
  return a;
}
/*
 * Get the modulus of a complex
 */
float floatComplexAbs(floatComplex x)
{
	float r;
	r=sqrt(pow(x.re,2)+pow(x.im,2));
	return r;
}
/* ADDITION FOR COMPLEX
 * parameters:
 * in)x--input complex x
 * y--intput complex y
 * out)z--output complex z
 */

floatComplex floatCAdd(floatComplex x,floatComplex y)
{
  floatComplex z;
  z.re=x.re+y.re;
  z.im=x.im+y.im;
  return z;
}
/*
 * MULTIPLICATION FOR COMPLEX
 * parameters:
 * in)x--input complex x
 * y--intput complex y
 * out)z--output complex z
 */

floatComplex floatCMpy(floatComplex x,floatComplex y)
{
  floatComplex z;
  z.re=x.re*y.re-x.im*y.im;
  z.im=x.re*y.im+x.im*y.re;
  return z;
}
/*
 *Multiplication for complex matrix
 *in)two matrixs:m*p,p*n
 *out)matrix:m*n
 */
floatComplexMatrix *floatComplexMatrixMpy(floatComplexMatrix *mat1,floatComplexMatrix *mat2)
{
	floatComplexMatrix *r;
	floatComplex temp;
	int i,j,k;
	int row,col;
	int n;
	
	row=mat1->row;
	col=mat2->col;
	n=mat1->col;
	r=floatComplexMatrixAlloc(row,col);
	
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
		{
			temp=floatComplexSet(0,0);
			for(k=0;k<n;k++)
				temp=floatCAdd(temp,floatCMpy(mat1->a[i][k],mat2->a[k][j]));
			r->a[i][j]=floatComplexSet(temp.re,temp.im);
		}
	return r;
}
/*
 * Get the Inverse matrix of a matrix
 */
floatComplexMatrix *floatComplexMatrixInverse(floatComplexMatrix *mat)
{
	floatComplexMatrix *p,*pTemp;
	floatComplexMatrix *q,*qTemp;
	floatComplexMatrix *e;
	floatComplex *rowTemp;
	int i,j,k,l;
	int n;
		
	n=mat->col;
	p=floatComplexMatrixAlloc(n,n);/* used for calculating the temp matrix after several kinds of preliminary transform */
	q=floatComplexMatrixAlloc(n,n);/* used for calculating the inverse */
	e=floatComplexMatrixAlloc(n,n);/* preliminary matrix */
	rowTemp=malloc(n*sizeof(floatComplex));
	/* Initializing p,q,e */
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			p->a[i][j]=mat->a[i][j];
			if(i==j)
			{
				q->a[i][j]=floatComplexSet(1,0);
				e->a[i][j]=floatComplexSet(1,0);
			}
			else
			{
				q->a[i][j]=floatComplexSet(0,0);
				e->a[i][j]=floatComplexSet(0,0);
			}
		}
	/*Assert the diagonal element is non-zero,o.w. exchange the row with another*/
	for(k=0;k<n;k++)
	{	
		if(floatComplexAbs(p->a[k][k])==0||floatComplexAbs(p->a[k][k])<1E-3)
		{
			for(l=k+1;l<n;l++)
			{
				if(floatComplexAbs(p->a[l][k])!=0&&floatComplexAbs(p->a[l][k])>1E-3)
				{
					for(i=0;i<n;i++)
					{
						rowTemp[i]=floatComplexSet(e->a[k][i].re,e->a[k][i].im);
						e->a[k][i]=floatComplexSet(e->a[l][i].re,e->a[l][i].im);
						e->a[l][i]=floatComplexSet(rowTemp[i].re,rowTemp[i].im);
					}
					break;
				}
			}
		}
		pTemp=floatComplexMatrixMpy(e,p);
		qTemp=floatComplexMatrixMpy(e,q);
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
			{
				p->a[i][j]=floatComplexSet(pTemp->a[i][j].re,pTemp->a[i][j].im);
				q->a[i][j]=floatComplexSet(qTemp->a[i][j].re,qTemp->a[i][j].im);
			}
		floatComplexMatrixFree(pTemp);
		floatComplexMatrixFree(qTemp);
		/*The diagonal element must be normalized to one by using preliminary row transform */
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
			{
				if(i==j)
				{
					e->a[i][j]=floatComplexSet(1,0);
				}
				else
				{
					e->a[i][j]=floatComplexSet(0,0);
				}
			}
		e->a[k][k]=floatComplexSet(p->a[k][k].re/pow(floatComplexAbs(p->a[k][k]),2),
					-p->a[k][k].im/pow(floatComplexAbs(p->a[k][k]),2));
		pTemp=floatComplexMatrixMpy(e,p);
		qTemp=floatComplexMatrixMpy(e,q);
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
			{
				p->a[i][j]=floatComplexSet(pTemp->a[i][j].re,pTemp->a[i][j].im);
				q->a[i][j]=floatComplexSet(qTemp->a[i][j].re,qTemp->a[i][j].im);
			}
		floatComplexMatrixFree(pTemp);
		floatComplexMatrixFree(qTemp);
		e->a[k][k]=floatComplexSet(1,0);
		/*The col where the k-th diagonal element exist must be transformed to 0 by using preliminary row transform */
		for(l=0;l<n;l++)
		{
			if(l!=k)
			{
				e->a[l][k]=floatComplexSet(-p->a[l][k].re,-p->a[l][k].im);
				pTemp=floatComplexMatrixMpy(e,p);
				qTemp=floatComplexMatrixMpy(e,q);
				for(i=0;i<n;i++)
					for(j=0;j<n;j++)
					{
						p->a[i][j]=floatComplexSet(pTemp->a[i][j].re,pTemp->a[i][j].im);
						q->a[i][j]=floatComplexSet(qTemp->a[i][j].re,qTemp->a[i][j].im);
					}
				floatComplexMatrixFree(pTemp);
				floatComplexMatrixFree(qTemp);
				e->a[l][k]=floatComplexSet(0,0);
			}
		}
	}
	floatComplexMatrixFree(p);
	floatComplexMatrixFree(e);
	free(rowTemp);
	return q;
}