线性方程组求解C语言实现

#include<iostream> #include<fstream> #include<math.h> #include"Gaus_Elimination.h" using namespace std; #define ROW 4 #define COL 4 //fabs(double a)//浮点型的绝对值是fabs,abs()只能用于整形 /*double ** New_GetMemory(int row,int col) { int i; double ** A; A = new double *[row+1]; for (i = 0;i <= row; ++ i) { A[i] = new double[col + 1]; memset(A[i],0,(col+1)*sizeof(double)); } return A; }*/ /* bool Gaus_Elim_with_Partial_Pivoting(double **Input,double **Target,double **Output,int row,int col) { int i,j,k; double ** A = New_GetMemory(row,col); double ** B = New_GetMemory(row,1); //检查内存分配的情况 if(A == NULL || B == NULL) return false; for (i = 1; i <= row ; ++ i)//保存传递过来的数据 { B[i][1] = Target[i][1]; for (j = 1; j <= col ; ++ j) { A[i][j] = Input[i][j]; } } //行变换 for (i = 1; i < row ; ++ i) { double temp ; if(fabs(A[i][i]) > 1e-7)//相对而言对角线上的数应该比较大 { for (j = i + 1; j <= row ; ++ j) { temp = A[j][i]/A[i][i]; for (k = i; k <= col ; ++ k) { A[j][k] -= temp * A[i][k]; ;//便于调试时加断点 } B[j][1] = B[j][1] - temp * B[i][1]; ; } } else//对角线上的元素为0或者很接近 { //Partial_Pivoting主要是解决对角线上元素为零的情况 //实质就是将为零或者接近于零的元素与该列中元素最大的值所在的行互换 int flag = i;//用于存储最大的值所在的行 double Max = fabs(A[i][i]);//用于保存最大值 for (j = i + 1; j <= row ; ++ j) { if (fabs(A[j][i])-Max > 1e-7) { Max = fabs(A[j][i]); flag = j; } } if(flag == i)//下面的值全部是零 return false; //实现互换 //系数矩阵的相互互换 for (k = 1 ; k <= col ; ++ k) { A[i][k] = A[i][k] + A[flag][k]; A[flag][k] = A[i][k] - A[flag][k]; A[i][k] = A[i][k] - A[flag][k]; } //矩阵b的互换 B[i][1] = B[i][1] + B[flag][1]; B[flag][1] = B[i][1] - B[flag][1]; B[i][1] = B[i][1] - B[flag][1]; i -= 1;//重新进行i行的变换 } } //反向消元，求各个元素的和 for(i = row; i >= 1; --i) { if(row == i) { if(fabs(A[i][i])>1e-15)//不等于0 Output[i][1] = B[i][1]/A[i][i]; else return false; } else { double sum = 0; for (j = i + 1; j <= col ; ++ j) //xii*aii+xi(i+1)*ai(i+1)+...+xi(i+k)ai(i+k) = b[i][1]; //xii = (b[i][1]-(xi(i+1)*ai(i+1)+...+xi(i+k)ai(i+k)))/aii; sum += A[i][j]*Output[j][1]; if(fabs(A[i][i])>1e-7) Output[i][1] = (B[i][1]-sum)/A[i][i]; else return false; } } for (i = 1; i <= row ; ++ i) { delete [] A[i]; delete [] B[i]; } delete [] A; delete [] B; A = NULL; B = NULL; return true; }*/ /* bool Gaus_Elim_with_Partial_Pivoting(double **Input,double **Target,double **Output,int row,int col) { int i,j,k; double ** A = New_GetMemory(row,col); double ** B = New_GetMemory(row,1); //检查内存分配的情况 if(A == NULL || B == NULL) return false; for (i = 1; i <= row ; ++ i)//保存传递过来的数据 { B[i][1] = Target[i][1]; for (j = 1; j <= col ; ++ j) { A[i][j] = Input[i][j]; } } //行变换 for (i = 1; i < row ; ++ i) { double temp ; //Partial_Pivoting主要是解决对角线上元素为零的情况 //实质就是将该列中最大值所在的行作为基准行，用于加减操作 int flag = i;//用于存储最大的值所在的行 double Max = fabs(A[i][i]);//用于保存最大值，fabs()主要是针对浮点型数据 for (j = i + 1; j <= row ; ++ j)//找到最大的主元 { if (fabs(A[j][i])-Max > 1e-7) { Max = fabs(A[j][i]); flag = j; } } //实现互换 //系数矩阵的相互互换 if(i != flag) //说明要交换。 //如果是i == flag，说明是同一组数，不能交换，因为每一次操作都是针对同一个地址 { for (k = 1 ; k <= col ; ++ k) { //典型的不引入变量参量值相互转换算法 A[i][k] = A[i][k] + A[flag][k]; A[flag][k] = A[i][k] - A[flag][k]; A[i][k] = A[i][k] - A[flag][k]; } //矩阵b的互换 B[i][1] = B[i][1] + B[flag][1]; B[flag][1] = B[i][1] - B[flag][1]; B[i][1] = B[i][1] - B[flag][1]; } if(fabs(A[i][i]) > 1e-7)//判断主元最大值大于0 { for (j = i + 1; j <= row ; ++ j) { temp = A[j][i]/A[i][i]; for (k = i; k <= col ; ++ k) { A[j][k] -= temp * A[i][k]; ;//便于调试时加断点 } B[j][1] = B[j][1] - temp * B[i][1]; ; } } // else // return false; } //反向消元，求各个元素的和 for(i = row; i >= 1; --i) { if(row == i) { if(fabs(A[i][i])>1e-7)//不等于0 Output[i][1] = B[i][1]/A[i][i]; else return false; } else { double sum = 0; for (j = i + 1; j <= col ; ++ j) //xii*aii+xi(i+1)*ai(i+1)+...+xi(i+k)ai(i+k) = b[i][1]; //xii = (b[i][1]-(xi(i+1)*ai(i+1)+...+xi(i+k)ai(i+k)))/aii; sum += A[i][j]*Output[j][1]; if(fabs(A[i][i])>1e-7) Output[i][1] = (B[i][1]-sum)/A[i][i]; else return false; } } for (i = 1; i <= row ; ++ i) { delete [] A[i]; delete [] B[i]; } delete [] A; delete [] B; A = NULL; B = NULL; return true; } */ /* void Read_Data_From_File(const char* name,double **data,int nrow,int ncol) { int row,col; double ** A; A = new double *[nrow+1]; ifstream InFile(name);//文件读数据，保存到动态内存中。 for (row = 1;row <= nrow ; ++ row) { A[row]=new double[ncol+1]; for (col = 1; col <= ncol ; ++ col) InFile >> A[row][col]; } InFile.close(); for (row = 1; row <= nrow ; ++ row)//将数据保存到新的数据空间 for (col = 1; col <= ncol ; ++ col) { data[row][col] = A[row][col]; ;//用于调试，便于观察数据 } for (row = 1;row <= nrow ; ++ row) delete[] A[row]; delete[] A; A = NULL; } */ int main() { int row ,col; double ** data = New_GetMemory(ROW,COL); double ** target = New_GetMemory(ROW,1); double ** answer = New_GetMemory(ROW,1); if(data == NULL) { cout<<"The Memory is not enough !!!"<<endl; return 1; } if (target == NULL) { cout<<"The Memory is not enough !!!"<<endl; return 1; } if (answer == NULL) { cout<<"The Memory is not enough !!!"<<endl; return 1; } //双斜杠主要是为了解决'\'为特殊字符的问题，'\\'为转义字符\; Read_Data_From_File("E:\\程序库\\基本算法C++\\Gaus_Elim_with_Partial_Pivot\\My_data_to_Process.dat",data,ROW,COL); Read_Data_From_File("E:\\程序库\\基本算法C++\\Gaus_Elim_with_Partial_Pivot\\My_data_to_Find.dat",target,ROW,1); for (row = 1;row <= ROW ; ++ row) { for (col = 1; col <= COL ; ++ col) cout<<data[row][col]<<"\t"; cout<<endl; } for (row = 1;row <= ROW ; ++ row) cout<<target[row][1]<<endl; Gaus_Elim_with_Partial_Pivoting(data,target,answer,ROW,COL); cout<< "The Answer of linear system is :"<<endl; for (row = 1; row <= ROW ; ++ row) cout<<answer[row][1]<<endl; for (row = 1; row <= ROW; ++ row) { delete [] data[row]; delete [] target[row]; delete [] answer[row]; } delete [] data; delete [] target; delete [] answer; data = NULL; target = NULL; answer = NULL; return 0; }