/*
算法没有做过什么改动，加上了内存释放的部分。
算法介绍
矩阵求逆在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;i for(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;j b[j]=c[j][p];
for(i=0;i for(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;i a[i][k-1]=a[i][k-1]/a[k-1][k-1];
for(i=k;i for(j=k;j a[i][j]=a[i][j]-a[i][k-1]*a[k-1][j];
for(i=k;i b[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;j sum=sum+a[i][j]*b[j];
b[i]=(b[i]-sum)/a[i][i];
sum=0;
}
for(j=0;j c[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;i for (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;i mid = mid * ( *(p+i*n+j%n) );
result+=mid;
}
for(m=0;m {
mid=1; //逆序相减
for(i=0,j=n-1-m+n; i mid=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;i for(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; j printf("%10G ",*p++);
printf("\n");
}
printf("--------------\n");
}
