/*

算法没有做过什么改动，加上了内存释放的部分。

算法介绍

矩阵求逆在3D程序中很常见，主要应用于求Billboard矩阵。按照定义的计算方法乘法运算，严重影响了性能。在需要大量Billboard矩阵运算时，矩阵求逆的优化能极大提高性能。这里要介绍的矩阵求逆算法称为全选主元高斯-约旦法。

高斯-约旦法(全选主元)求逆的步骤如下：

首先，对于 k 从 0 到 n - 1 作如下几步：

从第 k 行、第 k 列开始的右下角子阵中选取绝对值最大的元素，并记住次元素所在的行号和列号，在通过行交换和列交换将它交换到主元素位置上。

这一步称为全选主元。

m(k, k) = 1 / m(k, k)

m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k

m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k

m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k

最后，根据在全选主元过程中所记录的行、列交换的信息进行恢复，恢复的原则如下：在全选主元过程中，先交换的行(列)后进行恢复；

原来的行(列)交换用列(行)交换来恢复。

*/

#include   #include   #include   void   main(void)

{

int   n;

printf("请输入你所求解的逆矩阵的阶数：n=");

scanf("%d",&n);

float   choose_the_main(float   **a,float   *b,int   k,int   n);

void   input_output(float   **a,float   **c,float   **q,int   n);

void   course(float   **a,float   **c,float   **q,float   *b,int   n);

float   **a,*b,**c,**q;

int   i;

//没有free()，有内存泄露。

a=(float   **)malloc(sizeof(float   *)*n);

b=(float   *)malloc(sizeof(float   )*n);

c=(float   **)malloc(sizeof(float   *)*n);

q=(float   **)malloc(sizeof(float   *)*n);

for(i=0;i{

*(a+i)=(float   *)malloc(sizeof(float   )*n);

*(c+i)=(float   *)malloc(sizeof(float   )*n);

*(q+i)=(float   *)malloc(sizeof(float   )*n);

}

printf("您现在做的是%d阶矩阵的逆矩阵求解\n",n);

input_output(a,c,q,n);

course(a,c,q,b,n);

//释放内存

//ly-2009-10-31

for(i=0;i{

free(*(a+i));

*(a+i) = NULL;

free(*(c+i));

*(c+i) = NULL;

free(*(q+i));

*(q+i) = NULL;

}

free(*a);

free(b);

free(*c);

free(*q);

//

}

float   choose_the_main(float   **a,float   *b,int   k,int   n)

{

float   d,t;

int   l,i,j;

d=a[k-1][k-1];

l=k-1;

for(i=k;i{

if(fabs(a[i][k-1]>fabs(d)))

{

d=a[i][k-1];

l=i;

}

}

if(d==0)

printf("The   system   is   error!\n");

else

{

if(l!=k-1)

{

for(j=k-1;j{

t=a[l][j];

a[l][j]=a[k-1][j];

a[k-1][j]=t;

}

t=b[l];b[l]=b[k-1];b[k-1]=t;

}

}

return   (d);

}

/*a: 原始矩阵；c: 单位矩阵*/

void   input_output(float   **a,float   **c,float   **q,int   n)

{

int   i,j;

printf("请按行的顺序依次输入矩阵中元素的值(共%d项)：",n*n);

for(i=0;i<=n-1;i++)

for(j=0;j<=n-1;j++)

scanf("%f",&a[i][j]);

for(i=0;ifor(j=0;j<=n-1;j++)

{

if(i==j)   c[i][j]=1;

else   c[i][j]=0;

q[i][j]=a[i][j];

}

printf("您输入的(AE)矩阵为(利用(AE)-->(E(A的逆))：\n");

for(i=0;i<=n-1;i++)

{

for(j=0;j<=n-1;j++)

printf("%.5f     ",a[i][j]);

for(j=0;j<=n-1;j++)

printf("%.5f     ",c[i][j]);

printf("\n");

}

printf("\n");

}

void   course(float   **a,float   **c,float   **q,float   *b,int   n)

{

float   sum=0,h;

int   i,j,k,p,flag=0;

for(p=0;p{

for(j=0;jb[j]=c[j][p];

for(i=0;ifor(j=0;j<=n-1;j++)

{

a[i][j]=q[i][j];

}

for(k=1;k<=n-1;k++)

{

h=choose_the_main(a,b,k,n);  //选主元

if(h==0)

{

printf("因为矩阵的行列式的值为0，所以不能用此程序做！程序退出！\n");

flag=1;

break;

}

else

{

for(i=k;ia[i][k-1]=a[i][k-1]/a[k-1][k-1];

for(i=k;ifor(j=k;ja[i][j]=a[i][j]-a[i][k-1]*a[k-1][j];

for(i=k;ib[i]=b[i]-(a[i][k-1]*b[k-1]);

}

}

if(flag==1)   break;

if(h!=0)

{

b[n-1]=b[n-1]/a[n-1][n-1];

for(i=n-2;i>=0;i--)

{

for(j=i+1;jsum=sum+a[i][j]*b[j];

b[i]=(b[i]-sum)/a[i][i];

sum=0;

}

for(j=0;jc[j][p]=b[j];

}

}

if(flag==0)

{

printf("所求矩阵的逆矩阵为：\n");

for(i=0;i<=n-1;i++)

{

for(j=0;j<=n-1;j++)

printf("%.5f     ",c[i][j]);

printf("\n");

}

}

}

//8888888888888888888888888888888888888888888888888888888888888888888

//传统方法：

//***************************

//求任何一个矩阵的逆矩阵

//ly-09-10-27:求逆的传统方法，通过求代数余子式

//***************************

#include #include

void main( void )

{

float *buffer,*p;  //定义数组首地址指针变量

short int row,num; //定义矩阵行数row及矩阵元素个数

short int i,j;

float determ;     //定义矩阵的行列式

float comput_D(float *p,short int n);     //求矩阵的行列式

float Creat_M(float *p, short int m,short int n,short int k); //求代数余子式

void Print( float *p,short int n);    //打印n×n的矩阵

printf("\nPlease input the number of rows: ");

scanf("%d",&row);

num=2 * row * row;

buffer = (float *)calloc(num, sizeof(float));    //分配内存单元

p=buffer;

if(p != NULL)

{

for(i=0;i{

printf("Input the number of %d row ",i+1);

for(j=0;j{

scanf("%f",p++);

}

}

}

else

printf( "Can't allocate memory\n" );

printf("\nThe original matrix is:\n");

Print(buffer,row);    //打印该矩阵

determ=comput_D(buffer,row);    //求整个矩阵的行列式

p=buffer + row * row;

if (determ != 0)

{

for (i=0;ifor (j=0; j*(p+j*row+i)=  Creat_M(buffer,i,j,row)/determ;

printf("The determinant is %G\n",determ);

p=buffer + row * row;

printf("\nThe inverse matrix is:\n");

Print(p,row);    //打印该矩阵

}

else

printf("The determnant is 0, and there is no inverse matrix !\n");

free( buffer );

}

//--------------------------------------------------------

//功能：求矩阵 n X n 的行列式

//入口参数：矩阵首地址 p；矩阵行数 n

//返回值：矩阵的行列式值

//--------------------------------------------------------

float comput_D(float *p,short int n)

{

short int i,j,m;        //i--row; j--column

short int lop=0;

float result=0;

float mid=1;

if (n!=1)

{

lop=(n==2)?1:n;    //控制求和循环次数，若为2阶，则循环1次，否则为n次

for(m=0;m{

mid=1;         //顺序求和

for(i=0,j=m;imid = mid * ( *(p+i*n+j%n) );

result+=mid;

}

for(m=0;m{

mid=1;         //逆序相减

for(i=0,j=n-1-m+n; imid=mid * ( *(p+i*n+j%n));

result-=mid;

}

}

else result=*p;

return(result);

}

//----------------------------------------------------

//功能：求k×k矩阵中元素A(mn)的代数余子式

//入口参数：k×k矩阵首地址；元素A的下标m,n; 矩阵行数 k

//返回值： k×k矩阵中元素A(mn)的代数余子式

//----------------------------------------------------

float Creat_M(float *p, short int m,short int n,short int k)

{

short int len;

short int i,j;

float mid_result=0;

short int quo=1;

float *p_creat,*p_mid;

len=(k-1)*(k-1);

p_creat = (float *)calloc(len, sizeof(float));    //分配内存单元

p_mid=p_creat;

for(i=0;ifor(j=0;j{

if (i!=m && j!=n)

*p_mid++ =* (p+i*k+j);

}

//    Print(p_creat,k-1);

quo = (m + n) %2==0 ? 1:-1;

mid_result = (float ) quo * comput_D(p_creat,k-1);

free(p_creat);

return(mid_result);

}

//-------------------------------------------

//功能：打印n×n的矩阵

//入口参数：n×n矩阵的首地址；该矩阵的行数 n

//返回值： 无

//-------------------------------------------

void Print( float *p,short int n)

{

int i,j;

for (i=0;i{

for (j=0; jprintf("%10G ",*p++);

printf("\n");

}

printf("--------------\n");

}