2015年1月31日
高斯消元法求逆矩阵_C语言二维数组
高斯消元法是线性代数中的一个算法，可用来为线性方程组求解，以及求出可逆方阵的逆矩阵。在高斯消元法求逆矩阵C代码一文中讲述的逆矩阵算法用的是一维数组，本文将在C语言中用二维数组来实现。为什么要用二维数组来实现呢？做MCU的人都知道，在MCU中运行的程序是以微秒级来算，而且必须是争毫秒夺微秒。以100M主频的MCU而言，一维数组实现的高斯消元法执行一次需要40us，时间过长，严重影响程序的执行效率，所以必须要缩短这段程序的执行时间。
将高斯消元法求逆矩阵C代码的实现方式修改为二维数组后，执行时间(100M的MCU)可由原来的40us缩短为22us。对比可知，二维数组可以省去大量数组下标计算的时间。本文在Visual studio2008建立VC++ Win32控制台应用程序进行测试，测试结果如下：
将高斯消元法改为二维数组后，基本已经是该算法的最快速度，如果待求解的二维数组的阶数确定，可以将二维数组中的循环展开，大概可以再缩短2us左右(100M的MCU)。
|
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 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 |
#include "stdafx.h" #include <iostream> #include <stdlib.h> #include <math.h> #include <stdio.h> #include <time.h> #include<windows.h> using namespace std; unsigned char matrixInversion(); static double s_E[4][4]={{0.2368,0.2471,0.2568,1.2671}, {1.1161,0.1254,0.1397,0.1490}, {0.1582,1.1675,0.1768,0.1871}, {0.1968,0.2071,1.2168,0.2271}}; double eslaps; signed char is[4],js[4]; int _tmain(int argc, _TCHAR* argv[]) { int i,j,k; static double **b = new double *[4]; // 拷贝a static double unitMatrix[4][4]; for (i=0; i< 4; i++) { b[i] = new double[4]; for (j=0; j< 4; j++) b[i][j]=s_E[i][j]; // 拷贝a } LARGE_INTEGER nFreq; LARGE_INTEGER start,end; QueryPerformanceFrequency(&nFreq);//返回每秒嘀哒声的个数,即频率 QueryPerformanceCounter(&start); //获取开始时计数器的数值 unsigned char bOK = matrixInversion(); // 计算逆矩阵，结果在b中 QueryPerformanceCounter(&end); //获取结束时计数器的数值 eslaps=(double)(end.QuadPart-start.QuadPart)/(double)nFreq.QuadPart; //计数器滴答的次数与其频率的比值即为流逝的时间值 for (i = 0; i < 4; i ++) { for (j = 0; j < 4; j ++) { for(k = 0; k < 4; k ++) { unitMatrix[i][j] = unitMatrix[i][j] + s_E[i][k] * b[k][j]; } } } static int IntMaxtri[4][4]; for (int i = 0; i < 4; i ++) { for (int j = 0; j < 4; j ++) { IntMaxtri[i][j] = (int)(b[i][j] * 8192); } } if (i!=0) { printf("MAT A IS:\n"); for (i=0; i<=3; i++) { for (j=0; j<=3; j++) printf("%13.4f",s_E[i][j]); printf("\n"); } printf("\nMAT A- IS:\n"); for (i=0; i<=3; i++) { for (j=0; j<=3; j++) printf("%13.4f",b[i][j]); printf("\n"); } printf("\nMAT c- IS:\n"); for (i=0; i<=3; i++) { for (j=0; j<=3; j++) printf("%13.4f",unitMatrix[i][j]); printf("\n"); } printf("eslaps=%13.8f s\n",eslaps); } return 0; } unsigned char matrixInversion() { signed char i,j,k; double d,p; for(k = 0; k < 4; k ++) { is[k] = 0; js[k] = 0; d = 0.0; for (i = k; i <= 3; i ++) { for (j = k; j <= 3; j ++) { //将待求逆矩阵的数据放在二维数组s_E中，求得的逆矩阵也将存储在该数组中 p = fabs(s_E[i][j]); if(p>d) { d = p; is[k] = i; js[k] = j; } } } if ( 0.0 == d ) { return 0; } if (is[k] != k) { for(j = 0; j <= 3; j ++) { p = s_E[k][j]; s_E[k][j] = s_E[is[k]][j]; s_E[is[k]][j] = p; } } if (js[k] != k) { for (i = 0; i <= 3; i ++) { p = s_E[i][k]; s_E[i][k] = s_E[i][js[k]]; s_E[i][js[k]] = p; } } s_E[k][k] = 1.0/s_E[k][k]; for(j = 0; j <= 3; j ++) { if(j != k) { s_E[k][j] *= s_E[k][k]; } } for (i = 0; i <= 3; i ++) { if (i != k) { for(j = 0; j <= 3; j ++) { if(j != k) { s_E[i][j] -= s_E[i][k]*s_E[k][j]; } } } } for(i = 0; i <= 3; i ++) { if(i != k) { s_E[i][k] = -s_E[i][k]*s_E[k][k]; } } } for(k = 3; k >= 0; k --) { if(js[k] != k) { for(j = 0; j <= 3; j ++) { p = s_E[k][j]; s_E[k][j] = s_E[js[k]][j]; s_E[js[k]][j] = p; } if (is[k] != k) { for (i = 0; i <= 3; i ++) { p = s_E[i][k]; s_E[i][k] = s_E[i][is[k]]; s_E[i][is[k]] = p; } } } } return 1; } |
