利用初等行变换求矩阵的逆运算的代码.docVIP

/**************reverse matrix************************/ //利用的是AX=B，X=A’B,这里B=E；进行初等行变换求解，把左边化为单位阵，右边就是A矩阵的逆矩阵； void swap(double *a,int i,int line,int n) // exchange line//交换行位置，i控制行号，line也是行号，//n是矩阵列数 { int j; double temp; for(j=0;jn;j++) { temp=a[i*n+j]; a[i*n+j]=a[line*n+j]; a[line*n+j]=temp; } } void calculate(double *p,double *q,int n) //形成上三角阵 { int i,j,k,m,line; double max,temp,mmul; for(i=0;in;i++) { max=fabs(p[i*n+i]); temp=p[i*n+i]; line=i; for(j=i+1;jn;j++) //选出每一列中最大值，并且用line记录行号, //用temp记录该最大值，用j来控制行号 { if(fabs(p[j*n+i])max) { max=fabs(p[j*n+i]); temp=p[j*n+i]; line=j; } } if(max=1e-5) { printf(no inverse array\n); return; } if(line!=i) { swap(p,i,line,n);//将每一列中最大行换到i行 swap(q,i,line,n); } for(k=0;kn;k++) { p[i*n+k]/=temp;//将i行的每个数都除以i行i列的值，将i行i列化为1 q[i*n+k]/=temp; } for(k=i+1;kn;k++)//将i列i行下面其它行的值都化为0 { mmul=p[k*n+i]; for(m=0;mn;m++) { p[k*n+m]-=p[i*n+m]*mmul;//每一行都减去Ak0*i行对应列的值 q[k*n+m]-=q[i*n+m]*mmul; } } } } void backcalculate(double*p,double*q,int n)//形成单位矩阵 { int i,j,k; double mmul; for(i=n-1;i0;i--) { for(j=i-1;j=0;j--)//从下往上每一行进行计算，与前面相反 { mmul=p[j*n+i]; p[j*n+i]-=p[i*n+i]*mmul; for(k=0;kn;k++) { q[j*n+k]-=q[i*n+k]*mmul; } } } } void unitmatrix(double *q,int n) //形成单位矩阵，矩阵q是一个单位阵 { int i,j; for(i=0;in;i++) { for(j=0;jn;j++) { if(i==j) { q[i*n+j]=1; } else q[i*n+j]=0; } } } void MatrixInverse(double *A,double *C,int n)//矩阵求逆，A为要求的矩阵，C为求得的逆矩阵 { unitmatrix(C,n);//将C初始化为单位阵 calculate(A,C,n);//将A化为了上三角阵 backcalculate(A,C,n); //将A化为了单位阵，这时C即为所求 }