linear::CLinearEqualtion l(3);l.init_A("A.txt");l.init_b("b.txt");l.lu_decomposition();l.lup_solve();l.show_y();l.show_x();l.save_solution();

#include <iostream>
#include <vector>
#include <cassert>
#include <fstream>using namespace std;namespace linear
{
#define array_1(__type) std::vector<__type>
#define array_2(__type) std::vector<array_1(__type)>class CLinearEqualtion{/*使用方法:1. 定义计算方程组的类对象,并初始化其系数矩阵的大小linear::CLinearEqualtion l(3);2. 读取系数阵 Al.init_A("A.txt");3. 读取 bl.init_b("b.txt");4. 对系数矩阵进行lu分解l.lu_decomposition();5. 求解方程组l.lup_solve();6. 显示反向替换的解l.show_y();7. 显示正向替换的解l.show_x();8. 保存方程的解l.save_solution();A.txt1    5    42    0    35    8    2b.txt1295x[0] = -0.157895x[1] = -0.0526316x[2] = 3.10526*/private:array_2(double) A;array_2(double) L;array_2(double) U;array_1(double) b;array_1(int) p;array_1(double) x;array_1(double) y;public:CLinearEqualtion(int n){assert(n>1);A = array_2(double)(n);L = array_2(double)(n);U = array_2(double)(n);for (int i= 0; i < n;i++){A[i] = array_1(double)(n);L[i] = array_1(double)(n);U[i] = array_1(double)(n);}b = array_1(double)(n);x = array_1(double)(n);y = array_1(double)(n);p = array_1(int)(n);} // !CLinearEqualtion(int n)virtual ~CLinearEqualtion(){} // !virtual ~CLinearEqualtion()public:void init_A(array_2(double) _A){for (int i = 0; i < _A.size(); i++)A[i].assign(_A[i].begin(), _A[i].end());} // !void init_A(array_2(double) _A)void init_A(std::string _fileName){std::ifstream in(_fileName.c_str());int i = 0;int j = 0;while(!in.eof()){double e = 0;in>>e;std::cout<<e<<std::endl;A[i][j] = e;if (j==A.size() - 1){i++;j = 0;}else{j++;}}} // !void init_A(std::string _fileName)void init_b(std::string _fileName){std::ifstream in(_fileName.c_str());int i = 0;while(!in.eof()){double e = 0;in>>e;std::cout<<e<<std::endl;b[i++] = e;}} // !void init_b(std::string _fileName)void lu_decomposition(){int n = A.size(); // get array sizefor (int i = 0; i < n; i++)p[i] = i;for (int k = 0; k < n; k++){int max = 0;int k_ = k;// get max e in col k.for (int i = k; i < n; i++){if (fabs(A[i][k]) > max){max = fabs(A[i][k]);k_ = i;}}// make sure not all is zero.if (max == 0)assert("singular matrix");// swap cur row,max row.if (k != k_){swap(p[k], p[k_]);for (int i = 0; i < n; i++)swap(A[k][i], A[k_][i]);                    }for (int i = k + 1; i < n; i++){// gen vA[i][k] /= A[k][k];for (int j = k + 1; j < n; j++){A[i][j] -= A[i][k] * A[k][j];}}}init_LU();} // !void lu_decomposition()void init_LU(){int n = A.size(); // get array sizefor (int i = 0; i < n; i++){for (int j = 0; j < n; j++){if (i > j){L[i][j] = A[i][j];} else{U[i][j] = A[i][j];if (i == j)L[i][i] = 1;}}}} // !void init_LU()void lup_solve(){int n = A.size();int i = 0, j = 0;for (i = 0; i < n; i++){x[i] = 0;y[i] = 0;}for (i = 0; i < n; i++){y[i] = b[p[i]];for (j = 0; j < i; j++)y[i] -= L[i][j] * y[j];}for (i = n - 1; i >= 0; i--){double delt = y[i];for (j = n - 1; j > i; j--)delt -= U[i][j] * x[j];x[i] = delt / U[i][j];}} // !void lup_solve()void show_y(){int n = A.size();std::cout << "###" << std::endl;for (int i = 0; i < n; i++){std::cout << "y[" << i << "]=" << y[i] << std::endl;}} // !void show_y()void show_x(){int n = A.size();std::cout << "###" << std::endl;for (int i = 0; i < n; i++)std::cout << "x[" << i << "]=" << x[i] << std::endl;} // !void show_x()void save_solution(){int n = A.size();ofstream out("result.txt", ios::out);out << "-------------------------------------" << std::endl;std::cout << "-------------------------------------" << std::endl;for (int i = 0; i < n; i++){                out << "x[" << i << "] = " << x[i]<< std::endl;std::cout << "x[" << i << "] = " << x[i]<< std::endl;}out << "-------------------------------------" << std::endl;std::cout << "-------------------------------------" << std::endl;out.close();}};
}