C语言求解线性方程组AX=b

    线性代数的一个核心问题就是线性方程组的求解问题，包括：解的存在性（是否有解），解的唯一性（有解，是无穷解还是唯一解）。通常将线性方程组写成AX=b的形式，其中，A为系数矩阵，X为未知数（向量的形式），b表示目标向量，举例如下：

    求解上述线性方程组，常采用高斯消元法，将AX=b写成增广矩阵的形式，采用初等行变化得到上三角矩阵，从而判断，方程组解的存在性以及是否为唯一解。

    根据高斯消元法，采用C语言编程求解线性方程组，程序功能为：要求用户输入系数矩阵的大小，接着输入系数矩阵和目标向量，然后给出用户输入的线性方程组的增广矩阵形式，若方程无解，则输出方程无解；若方程有解，输出行化简后的上三角矩阵形式，并说明是有唯一解还是无穷解；程序可以实现线性方程组的循环输入，直到用户输入q退出！

    程序源代码如下，程序中含有7个子函数，最为重要的是gauss_elimination()和exchange_row()这两个子函数，其功能实现了增广矩阵的行化简，这两个矩阵的参数均是指针数组，所以行交换其实是通过交换指针指向来实现的，不必交换整行整列，提高了函数的效率，另外exchange_row()是将含最大主元元素的行提到首行位置，这样可以检测主元位置是否为0，且能提高数值计算的精度。

/*======================================================

 # Author:      Joker@HIT    NolanRobot@163.com
 # Filetype:    C source code
 # Environment: Linux & Ubuntu 14.04
 # Tool:        Vim & Gcc
 # Date:        Mon Aug 22 2016
 # Descprition: Solve the AX=b by the Guass Elimination!

========================================================*/

#include<stdio.h>
#include<math.h>

#define ROW 100
#define COL 100

void get_coefficient(double [][COL], int, int);
void get_vector(double [], int);
void create_Ab(double [][COL], double [], int, int);
void show_matrix(double [][COL], int, int);
int guass_elimination(double *[ROW], int, int);
void exchange_row(double *[ROW], int, int);
void show_solution(double *[ROW], int, int);

int main(void)
{
	double Receptacle[ROW][COL];
	double Vector[ROW];
	double *Ab_pointer[ROW];
	int row, col, i, result;

	printf("Enter the coefficient matrix's size (less than %d * %d): ", ROW, COL-1);
	while(scanf("%d%d", &row, &col)==2){
		get_coefficient(Receptacle, row, col);
		get_vector(Vector, row);
		create_Ab(Receptacle, Vector, row, col);

		printf("The linear equations in the form of augmented matrix as follow:\n");
		show_matrix(Receptacle, row, col+1);
		
		for(i=0;i<ROW;i++)
			Ab_pointer[i] = Receptacle[i];

		result = guass_elimination(Ab_pointer, row, col+1);
		if(result==1){
			printf("There are infinite solutions for the linear equations!\n");
			show_solution(Ab_pointer, row, col+1);
		}
		else if(result==0){
			printf("There is only one solution for the linear equations!\n");
			show_solution(Ab_pointer, row, col+1);
		}
		else
			printf("There is no solution for the linear equations!\n");
			
		printf("Enter the coefficient matrix's size (less than %d * %d, q to quit): ", ROW, COL-1);
	}
	printf("Bye!\n");

	return 0;
}

void get_coefficient(double matrix[ROW][COL], int row, int col)
{
	int i, j;

	printf("Enter the coefficient (%d * %d):\n", row, col);
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
			scanf("%lf", &matrix[i][j]);

	return;
}

void get_vector(double vector[ROW], int n)
{
	int i;

	printf("Enter the vector (size is %d):\n", n);
	for(i=0;i<n;i++)
		scanf("%lf", &vector[i]);

	return;
}

void create_Ab(double matrix[ROW][COL], double vector[ROW], int row, int col)
{
	int i;

	for(i=0;i<row;i++)
		matrix[i][col] = vector[i];

	return;
}

void show_matrix(double matrix[ROW][COL], int row, int col)
{
	int i, j;

	for(i=0;i<row;i++){
		for(j=0;j<col;j++)
			printf("%-8.3f", matrix[i][j]);
		putchar('\n');
	}

	return;
}

int guass_elimination(double *matrix[ROW], int row, int col)
{
	int result, i, j, k;
	double coe;

	for(i=0;i<row-1;i++){
		exchange_row(matrix, i, row);
		if(fabs(*(matrix[i]+i))<0.00001)
			continue;
		for(j=i+1;j<row;j++){
			coe = *(matrix[j]+i) / *(matrix[i]+i);
			for(k=i;k<col;k++)
				*(matrix[j]+k) -= coe * *(matrix[i]+k);
		}
	}

	if(col-1>row)
		result = 1;
	else if(col-1==row){
		if(fabs(*(matrix[row-1]+row-1))>0.00001)
			result = 0;
		else{
			if(fabs(*(matrix[row-1]+row))>0.00001)
				result = -1;
			else
				result = 1;
		}
	}
	else{
		result = 0;
		for(i=0;i<row;i++)
			if(fabs(*(matrix[i]+col-2))<0.00001&&fabs(*(matrix[i]+col-1))>0.00001){
				result = -1;
				break;
			}
	}


	return result;
}

void exchange_row(double *matrix[ROW], int flag, int row)
{
	int i;
	double *temp;

	for(i=flag+1;i<row;i++)
		if(fabs(*(matrix[flag]+flag))<fabs(*(matrix[i]+flag))){
			temp = matrix[flag];
			matrix[flag] = matrix[i];
			matrix[i] = temp;
		}

	return;
}

void show_solution(double *matrix[ROW], int row, int col)
{
	int i, j;

	for(i=0;i<row;i++){
		for(j=0;j<col;j++)
			printf("%-8.3f", *(matrix[i]+j));
		putchar('\n');
	}

	return;
}

运行结果如下：

上述运行结果验证了3种情况，均正确，如无解直接给出结论，若有解（唯一解或无穷解），给出相应的上三角矩阵和相应的结论！

读者若对代码有什么建议，欢迎留言或邮箱探讨！