高斯消元法解向量方程Ax=b

#include <iostream>#include <malloc.h>#include <math.h>//非齐次线性方程组Ax=bdouble **mat_A; //存放系数矩阵Adouble *mat_b; //存放右值数组bdouble *result; //结果数组const int M=3; //方阵的阶const int N=3;///////////////////////////////////////void init_ary() { //初始化系数矩阵 A 和右值数组 b; std::cout << "系数矩阵的行为 " << M << "列为 " << N << '\n';std::cout << "input the mat_A and mat_b" << '\n';mat_A = (double **)malloc(M*sizeof(double));for(int i=0; i<M; i++)*(mat_A+i) = (double *)malloc(N*sizeof(double));std::cout << "input the mat_A " << '\n';for(i=0; i<M; i++)for(int j=0; j<N; j++)std::cin >> *(*(mat_A+i)+j);std::cout << "input the mat_b" << '\n';mat_b = (double *)malloc(M*sizeof(double));for(i=0; i<M; i++)std::cin >> *(mat_b+i);}/////////////////////////////////////////////void guassi() { //对增广矩阵进行高斯消元, 将系数矩阵变成上三角矩阵double *temp; //当该行主元为0时需要与主元列元素不为0的行交换, 临时数组temp用于交换temp = (double *)malloc(N*sizeof(double));/*----高斯消元游戏开始了----*/for(int i=0; i<M-1; i++) { //需要对前 M-1 行进行消元double keyData = *(*(mat_A+i)+i); //记录当前行的主元int keyRow = i; //记录当前主元行/*---------当前主元为0,需要找到新的主元行进行交换--------*/if(keyData == 0) { for(int k=i+1; k<M; k++)if(*(*(mat_A+k)+i) != 0) {keyRow = k;break; //找到当前主元列中当前元素不为0的行 k ,用于交换使其成成为真正的主元行}/*----交换两行-----------------------------*/for(int j=i; j<N; j++) {*(temp+j) = *(*(mat_A+i)+j);*(*(mat_A+i)+j) = *(*(mat_A+keyRow)+j);*(*(mat_A+keyRow)+j) = *(temp+j);}/*----交换两行-----------------------------*///交换后主元行元素更新, 故要重置keyRow,keyDatakeyRow = i; keyData = *(*(mat_A+i)+i); }/*---------当前主元为0,需要找到新的主元行进行交换--------*//*------开始对当前行进行高斯消元----------*///将主元列中主元以下的所有元素置0for(int k=i+1; k<M; k++) { //i 为主元行行号, 主元以下消元, 故从 i+1 开始double elm = *(*(mat_A+k)+i)/keyData; //主元实际就是对角线元素, 故主元所在列等于行值, 故主元列为 i; elm是用于消元的除式因子/*----本行所有元素进行消元相关变化---*/for(int j=i+1; j<N; j++) {*(*(mat_A+k)+j) -= *(*(mat_A+i)+j)*elm; //本行元素减去主元行元素乘以除式因子elm变为新的本行元素; i 是主元行的行号}*(mat_b+k) -= *(mat_b+i)*elm; //右值列向量矩阵也要发生相应变化 /*----本行所有元素进行消元相关变化---*/}/*------开始对当前行进行高斯消元----------*/}/*----高斯消元游戏结束了----*/}//////////////////////////////////////////////////////////////void cal_mat() { //求解Ax=bresult = (double *)malloc(N*sizeof(double));*(result+N-1) = (*(mat_b+M-1))/(*(*(mat_A+M-1)+N-1));int k=M-2; //标记mat_b的下标for(int i=N-2; i>=0; i--) { //回代求解double sum = 0.0;for(int j=N-1; j>i; j--) sum += *(*(mat_A+i)+j)* (*(result+j));*(result+i) = (*(mat_b+k)-sum)/(*(*(mat_A+i)+i));  k--;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            }std::cout << "out put the result\n";for(i=0; i<N; i++)std::cout << *(result+i) << "  ";std::cout << '\n';}//////////////////////////////////////////////////////////void main() {init_ary();guassi();cal_mat();}

 验证: