C语言求逆矩阵（用c语言求逆矩阵）

一般求逆矩阵的方法有两种，伴随阵法和初等变换法。但是这两种方法都不太适合编程。伴随阵法的计算量大，初等变换法又难以编程实现。
适合编程的求逆矩阵的方法如下：

对可逆矩阵A进行QR分解：A=QR
求上三角矩阵R的逆矩阵
求出A的逆矩阵：A^(-1)=R^(-1)Q^(H)

以上三步都有具体的公式与之对应，适合编程实现。
C语言实现代码：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 #include <stdio.h> #include <math.h> #define SIZE 8 double b[SIZE][SIZE]={0};//应该读作“贝尔塔”，注释中用B表示 double t[SIZE][SIZE]={0};//求和的那项 double Q[SIZE][SIZE]={0};//正交矩阵 double QH[SIZE][SIZE]={0};//正交矩阵的转置共轭 double R[SIZE][SIZE]={0};// double invR[SIZE][SIZE]={0};//R的逆矩阵 double invA[SIZE][SIZE]={0};//A的逆矩阵，最终的结果 //={0};// double matrixR1[SIZE][SIZE]={0}; double matrixR2[SIZE][SIZE]={0}; //double init[3][3]={3,14,9,6,43,3,6,22,15}; double init[8][8]={ 0.0938 , 0.5201 , 0.4424 , 0.0196 , 0.3912 , 0.9493 , 0.9899 , 0.8256, 0.5254 , 0.3477 , 0.6878 , 0.3309 , 0.7691 , 0.3276 , 0.5144 , 0.7900, 0.5303 , 0.1500 , 0.3592 , 0.4243 , 0.3968 , 0.6713 , 0.8843 , 0.3185, 0.8611 , 0.5861 , 0.7363 , 0.2703 , 0.8085 , 0.4386 , 0.5880 , 0.5341, 0.4849 , 0.2621 , 0.3947 , 0.1971 , 0.7551 , 0.8335 , 0.1548 , 0.0900, 0.3935 , 0.0445 , 0.6834 , 0.8217 , 0.3774 , 0.7689 , 0.1999 , 0.1117, 0.6714 , 0.7549 , 0.7040 , 0.4299 , 0.2160 , 0.1673 , 0.4070 , 0.1363, 0.7413 , 0.2428 , 0.4423 , 0.8878 , 0.7904 , 0.8620 , 0.7487 , 0.6787 }; /*/ 函数名:int main() 输入: 输出: 功能:求矩阵的逆 pure C language 首先对矩阵进行QR分解之后求上三角矩阵R的逆阵最后A-1=QH*R-1,得到A的逆阵。 作者：HLdongdong *////////////////////////////////////////////////////////////////////// int main() { int i;//数组行 int j;//数组列 int k;//代表B的角标 int l;//数组列 double dev; double numb;//计算的中间变量 double numerator,denominator; double ratio; /////////////////求B///////////////// for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { b[j][i]=init[j][i]; } for(k=0;k<i;++k) { if(i) { numerator=0.0; denominator=0.0; for(l=0;l<SIZE;++l) { numerator+=init[l][i]*b[l][k]; denominator+=b[l][k]*b[l][k]; } dev=numerator/denominator; t[k][i]=dev; for(j=0;j<SIZE;++j) { b[j][i]-=t[k][i]*b[j][k];//t init =0 !!! } } } } ///////////////////对B单位化,得到正交矩阵Q矩阵//////////////////// for(i=0;i<SIZE;++i) { numb=0.0; for(j=0;j<SIZE;++j) { numb+=(b[j][i]*b[j][i]); } dev=sqrt(numb); for(j=0;j<SIZE;++j) { Q[j][i]=b[j][i]/dev; } matrixR1[i][i]=dev; } /////////////////////求上三角R阵/////////////////////// for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { if(j<i) { matrixR2[j][i]=t[j][i]; } else if(j==i) { matrixR2[j][i]=1; } else { matrixR2[j][i]=0; } } } mulMatrix(matrixR1,matrixR2,SIZE,SIZE,SIZE,R); ///////////////////////QR分解完毕////////////////////////// printf("QR分解：\n"); printf("Q=\n"); for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { printf("%2.4f ",Q[i][j]); // } printf("\n"); } printf("R=\n"); for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { printf("%2.4f ",R[i][j]); // } printf("\n"); } /////////////////////求R的逆阵////////////////////////// for(i=SIZE-1;i>=0;--i) { invR[i][i]=1/R[i][i]; //R[i][i]=invR[i][i]; if(i!=(SIZE-1))//向右 { for(j=i+1;j<SIZE;++j) { invR[i][j]=invR[i][j]*invR[i][i]; R[i][j]=R[i][j]*invR[i][i]; } } if(i)//向上 { for(j=i-1;j>=0;--j) { ratio=R[j][i]; for(k=i;k<SIZE;++k) { invR[j][k]-=ratio*invR[i][k]; R[j][k]-=ratio*R[i][k]; } } } } /////////////////////////////////////////////////////// printf("inv（R）=\n"); for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { printf(" %2.4f ",invR[i][j]); // } printf("\n"); } ////////////////////结果和MATLAB差一个负号，神马鬼？？？？？？？？///////////////////// /////////////////////求QH////////////////////////// for(i=0;i<SIZE;++i)//实矩阵就是转置 { for(j=0;j<SIZE;++j) { QH[i][j]=Q[j][i]; } } ///////////////////////求A的逆阵invA///////////////////////////// mulMatrix(invR,QH,SIZE,SIZE,SIZE,invA); printf("inv（A）=\n"); for(i=0;i<SIZE;++i) { for(j=0;j<SIZE;++j) { printf(" %2.4f ",invA[i][j]); // } printf("\n"); } ///////////////////////结果与MATLAB的结果在千分位后有出入，但是负号都是对的^v^/////////////////////////// return 0; }

另附上矩阵乘法的子函数

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 /*/ 函数名:void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE]) 输入:依次是左矩阵，右矩阵，左矩阵高度，左矩阵宽度，右矩阵宽度，输出矩阵 输出: 功能:矩阵乘法 作者：HLdongdong *// void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE]) { int i,j,k; for(i=0;i<high1;++i) { for(j=0;j<weight2;j++) { for(k=0;k<weight;++k) { mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j]; } } } }