根据高斯消元法的思路，编写求逆矩阵的程序。

矩阵求逆的快速算法

算法介绍

矩阵求逆在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

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

实现(4阶矩阵)

float Inverse(CLAYMATRIX& mOut, const CLAYMATRIX& rhs)

{

CLAYMATRIX m(rhs);

DWORD is[4];

DWORD js[4];

float fDet = 10f;

int f = 1;

for (int k = 0; k < 4; k ++)

{

// 第一步，全选主元

float fMax = 00f;

for (DWORD i = k; i < 4; i ++)

{

for (DWORD j = k; j < 4; j ++)

{

const float f = Abs(m(i, j));

if (f > fMax)

{

fMax = f;

is[k] = i;

js[k] = j;

}

}

}

if (Abs(fMax) < 00001f)

return 0;

if (is[k] != k)

{

f = -f;

swap(m(k, 0), m(is[k], 0));

swap(m(k, 1), m(is[k], 1));

swap(m(k, 2), m(is[k], 2));

swap(m(k, 3), m(is[k], 3));

}

if (js[k] != k)

{

f = -f;

swap(m(0, k), m(0, js[k]));

swap(m(1, k), m(1, js[k]));

swap(m(2, k), m(2, js[k]));

swap(m(3, k), m(3, js[k]));

}

// 计算行列值

fDet = m(k, k);

// 计算逆矩阵

// 第二步

m(k, k) = 10f / m(k, k);

// 第三步

for (DWORD j = 0; j < 4; j ++)

{

if (j != k)

m(k, j) = m(k, k);

}

// 第四步

for (DWORD i = 0; i < 4; i ++)

{

if (i != k)

{

for (j = 0; j < 4; j ++)

{

if (j != k)

m(i, j) = m(i, j) - m(i, k) m(k, j);

}

}

}

// 第五步

for (i = 0; i < 4; i ++)

{

if (i != k)

m(i, k) = -m(k, k);

}

}

for (k = 3; k >= 0; k --)

{

if (js[k] != k)

{

swap(m(k, 0), m(js[k], 0));

swap(m(k, 1), m(js[k], 1));

swap(m(k, 2), m(js[k], 2));

swap(m(k, 3), m(js[k], 3));

}

if (is[k] != k)

{

swap(m(0, k), m(0, is[k]));

swap(m(1, k), m(1, is[k]));

swap(m(2, k), m(2, is[k]));

swap(m(3, k), m(3, is[k]));

}

}

mOut = m;

return fDet f;

}

比较

原算法 原算法(经过高度优化) 新算法

加法次数 103 61 39

乘法次数 170 116 69

需要额外空间 16 sizeof(float) 34 sizeof(float) 25 sizeof(float)

结果不言而喻吧。

我以前写过求逆矩阵的程序。不过没有用到结构体，你看看如何。

#include<stdioh>

void main()

{

int N;

printf("输入不超过10的矩阵的阶数N:\n");

scanf("%d",&N);

float a[10][10],b[10][20],c[10][10],t;

int i,j,m;

printf("请输入行列式不为0的矩阵A(%d阶):\n",N); //矩阵A的各元素存入二维数组a中。

for(i=0;i<N;i++)

for(j=0;j<N;j++)

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

//增广矩阵（A|E）存入二维数组b中

for(i=0;i<N;i++)

for(j=0;j<N;j++)

b[i][j]=a[i][j];

for(i=0;i<N;i++)

for(j=N;j<2N;j++)

b[i][j]=0;

for(i=0;i<N;i++)

b[i][N+i]=1;

for(m=0;m<N;m++) //对每行进行处理。

{

t=b[m][m]; //预存b[m][m]。

i=m;

while(b[m][m]==0)

{

b[m][m]=b[i+1][m];

i++;

}

if(i>m)

{

b[i][m]=t; //实现交换。

//交换其它各列相应位置的元素

for(j=0;j<m;j++)

{

t=b[m][j];

b[m][j]=b[i][j];

b[i][j]=t;

}

for(j=m+1;j<2N;j++)

{

t=b[m][j];

b[m][j]=b[i][j];

b[i][j]=t;

}

}

for(i=m+1;i<N;i++)

for(j=2N-1;j>=m;j--)

b[i][j]-=b[i][m]b[m][j]/b[m][m]; //m=0时，将第一行的-b[i][0]/b[0][0]倍加到以下各行。这样以下每行第一个元素b[i][0]就为0。

for(j=2N-1;j>=m;j--)

b[m][j]/=b[m][m]; //对第m行作行变换，同除以b[m][m]，使b[m][m]为1。

}

printf("第一步变换后得到的增广矩阵为：\n");

for(i=0;i<N;i++)

{

for(j=0;j<2N;j++)

printf("%35f ",b[i][j]);

printf("\n"); //实现了：每个i对应一个换行。

}

m=N-1;

while(m>0)

{

for(i=0;i<m;i++)

for(j=2N-1;j>=m;j--) //千万注意，此处j必须递减，否则b[i][m]先变为0，后面的计算就无效！

b[i][j]-=b[i][m]b[m][j];

m--;

}

printf("最后得到的增广矩阵为：\n");

for(i=0;i<N;i++)

{

for(j=0;j<2N;j++)

printf("%35f ",b[i][j]);

printf("\n"); //实现了：每个i对应一个换行。

}

for(i=0;i<N;i++) //将逆矩阵存入二维数组c中。

for(j=0;j<N;j++)

c[i][j]=b[i][N+j];

printf("故逆矩阵为：\n");

for(i=0;i<N;i++)

{

for(j=0;j<N;j++)

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

printf("\n"); //实现了：每个i对应一个换行。

}

}

! aa为原矩阵，b为存放aa的逆矩阵，n为矩阵aa的维数

subroutine nizhen(aa,b,n)

integer n,i,j,k

real:: aa(n,n),b(n,n),a(n,n)

a=aa

do i=1,n

b(i,i)=1

enddo

do i=1,n

b(i,:)=b(i,:)/a(i,i)

a(i,i:n)=a(i,i:n)/a(i,i)

do j=i+1,n

do k=1,n

b(j,k)=b(j,k)-b(i,k)a(j,i)

enddo

a(j,i:n)=a(j,i:n)-a(i,i:n)a(j,i)

enddo

enddo

do i=n,1,-1

do j=i-1,1,-1

do k=1,n

b(j,k)=b(j,k)-b(i,k)a(j,i)

enddo

enddo

enddo

end

楼上的是矩阵的转置，不是求逆矩阵！

c语言的方阵求逆的一个完整的程序如下，在tc20下运行良好：

#include <stdioh>

#include <stdlibh>

#include <conioh>

#include <mathh>

#define MAX 255

void MatrixMul(a,b,m,n,k,c) /实矩阵相乘/

int m,n,k; /m:矩阵A的行数, n:矩阵B的行数, k:矩阵B的列数/

double a[],b[],c[]; /a为A矩阵, b为B矩阵, c为结果，即c = AB /

{

int i,j,l,u;

/逐行逐列计算乘积/

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

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

{

u=ik+j; c[u]=00;

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

c[u]=c[u]+a[in+l]b[lk+j];

}

return;

}

int brinv(a,n) /求矩阵的逆矩阵/

int n; /矩阵的阶数/

double a[]; /矩阵A/

{

int is,js,i,j,k,l,u,v;

double d,p;

is=malloc(nsizeof(int));

js=malloc(nsizeof(int));

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

{

d=00;

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

/全选主元，即选取绝对值最大的元素/

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

{

l=in+j; p=fabs(a[l]);

if (p>d) { d=p; is[k]=i; js[k]=j;}

}

/全部为0，此时为奇异矩阵/

if (d+10==10)

{

free(is); free(js); printf(" >> This is a singular matrix, can't be inversed!\n");

return(0);

}

/行交换/

if (is[k]!=k)

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

{

u=kn+j; v=is[k]n+j;

p=a[u]; a[u]=a[v]; a[v]=p;

}

/列交换/

if (js[k]!=k)

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

{

u=in+k; v=in+js[k];

p=a[u]; a[u]=a[v]; a[v]=p;

}

l=kn+k;

a[l]=10/a[l]; /求主元的倒数/

/ a[kj]a[kk] -> a[kj] /

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

if (j!=k)

{

u=kn+j; a[u]=a[u]a[l];

}

/ a[ij] - a[ik]a[kj] -> a[ij] /

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

if (i!=k)

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

if (j!=k)

{

u=in+j;

a[u]=a[u]-a[in+k]a[kn+j];

}

/ -a[ik]a[kk] -> a[ik] /

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

if (i!=k)

{

u=in+k; a[u]=-a[u]a[l];

}

}

for (k=n-1; k>=0; k--)

{

/恢复列/

if (js[k]!=k)

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

{

u=kn+j; v=js[k]n+j;

p=a[u]; a[u]=a[v]; a[v]=p;

}

/恢复行/

if (is[k]!=k)

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

{

u=in+k; v=in+is[k];

p=a[u]; a[u]=a[v]; a[v]=p;

}

}

free(is); free(js);

return(1);

}

print_matrix(a,n)/打印的方阵a的元素/

int n; /矩阵的阶数/

double a[]; /矩阵a/

{

int i,j;

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

{

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

printf("%137f\t",a[in+j]);

printf("\n");

}

}

main()

{

int i,j,n=0;

double A[MAX],B[MAX],C[MAX];

static double a[4][4]={ {02368,02471,02568,12671},

{11161,01254,01397,01490},

{01582,11675,01768,01871},

{01968,02071,12168,02271}};

static double b[4][4],c[4][4];

clrscr();

puts("");

puts(" This program is to inverse a square matrix A(nxn) ");

puts("");

while(n<=0)

{

printf(" >> Please input the order n of the matrix (n>0): ");

scanf("%d",&n);

}

printf(" >> Please input the elements of the matrix one by one:\n >> ");

for(i=0;i<nn;i++)

{

scanf("%lf",&A[i]);

B[i]=A[i];

}

for(i=0;i<4;i++)

for(j=0;j<4;j++)

b[i][j]=a[i][j];

i=brinv(A,n);

if (i!=0)

{

printf(" Matrix A:\n");

print_matrix(B,n);

printf("\n");

printf(" A's Inverse Matrix A-:\n");

print_matrix(A,n);

printf("\n");

printf(" Product of A and A- :\n");

MatrixMul(B,A,n,n,n,C);

print_matrix(C,n);

}

printf("\n Press any key to quit");

getch();

}

需要在MATLAB中输入矩阵A： A=[1 2 3;2 2 1;3 4 3]，回车；

输入：inv(A)或A^-1；回车。

>> syms a b c d (定义变量)

>> A=[a,b;c,d] (定义矩阵)

A =

[ a,b]

[ c,d]

>> inv(A) (求矩阵的逆)

ans = (结果)

[ d/(ad - bc),-b/(ad - bc)]

[ -c/(ad - bc),a/(ad - bc)]

MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

matalb 有计算逆矩阵的函数 inv()

例如：

>> A=magic(3)

A =

8 1 6

3 5 7

4 9 2

>> inv(A)

ans =

01472 -01444 00639

-00611 00222 01056

-00194 01889 -01028

以上就是关于根据高斯消元法的思路，编写求逆矩阵的程序。全部的内容，包括:根据高斯消元法的思路，编写求逆矩阵的程序。、C语言用二维数组实现矩阵求逆、fortran语言矩阵求逆等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！