高斯消元法解方程Ax=b的C源代码实现

采用列选主元方式实现，过程中维护了两个数组，分别记录改行是否已经选取主元，与选取主元行数的顺序。回代时利用选取顺序逆序得到最终结果。输入时先输入矩阵A的维数（方阵）n，再分别输入n个行向量；最后输入b。代码如下（由于使用了_Bool类型数组，需添加编译选项-std=c99，或者你可以直接改成int类型的：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
double Abs(double x) { return (x > 0) ? x : -x; }
int* GaussE(double **A, double *b, int size);
double* backSub(double **A, double *b, int *order, int size);
int findMaxLin(double **A, _Bool *flag, int col, int size);
int main(void) {
int size;
double **A;
double *b;
double *slt;
int *order;
int cnt;
scanf("%d", &size);
getchar();
A = (double **)malloc(sizeof(double *)*size);
for (cnt = 0; cnt < size; cnt++) {
A[cnt] = (double *)malloc(sizeof(double)*size);
}
b = (double *)malloc(sizeof(double)*size);

for (cnt = 0; cnt < size*size; cnt++) {
scanf("%lf", &A[cnt / size][cnt - cnt / size*size]);
}
for (cnt = 0; cnt < size; cnt++) {
scanf("%lf", b + cnt);
}

order = GaussE(A, b, size);
slt = backSub(A, b, order, size);

//print the results
printf ("Matrix after Gauss Elimination:\n");
for (cnt = 0; cnt < size*size; cnt++) {
printf("%.3lf ", A[cnt / size][cnt - cnt / size*size]);
if ((cnt + 1) % size == 0)
printf("\n");
}
printf ("\n-------------------------------------\n");
printf ("Solution = (");
for (cnt = 0; cnt < size-1; cnt++) {
printf("%.3lf, ", slt[cnt]);
}
printf ("%.3lf)\n", slt[cnt]);

//free the memory
for (cnt = 0; cnt < size; cnt++) {
free(A[cnt]);
}
free(A);
free(b); free(order); free(slt);
return 0;
}
int* GaussE(double **A, double *b, int size) {
int i, j, k;
int maxLinNum;
int *order = (int *)malloc(sizeof(int)*size);
double division;
_Bool *flag = (_Bool *)malloc(sizeof(_Bool)*size);  //Mark the line vector of which the first element has been used as the pivot
memset(flag, 0, sizeof(_Bool)*size);
for (i = 0; i < size; i++) {
order[i] = findMaxLin(A, flag, i, size);        //Record the order of elimination, for the convenience of back substitution
for (j = 0; j < size; j++) {                    //Go thru all the line vectors
if (flag[j] == 0) {
division = A[j][i] / A[order[i]][i];        //get the division: pivot/A[j][i]
for (k = i + 1; k < size; k++) {
A[j][k] -= A[order[i]][k] * division;
}
A[j][i] = 0;
b[j] -= b[order[i]] * division;
}
else
continue;
}
}
free(flag);
return order;
}
int findMaxLin(double **A, _Bool *flag, int col, int size) {
//Choose the most suitable pivot.
double max = 0;
int cnt;
int maxIndex;
for (cnt = 0; cnt < size; cnt++) {
if (max == 0 && flag[cnt] == 0) {
max = Abs(A[cnt][col]);
maxIndex = cnt;
continue;
}
else if (flag[cnt] == 0 && Abs(A[cnt][col]) > max) {
max = Abs(A[cnt][col]);
maxI
4000
ndex = cnt;
continue;
}
}
flag[maxIndex] = 1;
return maxIndex;
}
double* backSub(double **A, double *b, int *order, int size) {
int cnt;
int beg;
double *slt = (double *)malloc(sizeof(double)*size);
for (cnt = size-1; cnt >= 0; cnt--) {
for (beg = cnt+1; beg < size; beg++) {
b[order[cnt]] -= A[order[cnt]][beg]*slt[beg];
}
slt[cnt] = b[order[cnt]]/A[order[cnt]][cnt];
}
return slt;
}