解线性方程组.docVIP

//线性方程组列主元高斯消元法求解 #include stdafx.h #includememory #if _DEBUG #includeiostream //测试时使用 #endif using namespace std; //绝对值函数 __inline double _abs(double v) { return v 0 ? -v : v; } #if _DEBUG void test(int n, const double *pd) { int i, j, k = 0; for (i = 0; i n; i++) { for (j = 0; j = n; j++) { cout pd[k] ; k++; } cout endl; } } #endif //线性方程组列主元高斯消元法 //n 方程元数；pCoef 系数，必须以行主序方式存放的二维数组； //pOut 长度为 n 的一维数组(调用者负责维护)，用于输出数据 //返回值：0 成功，-1 无解，1 申请内存失败, 2 不定解。 int GaussJordanElimination(int n, const double *pCoef, double *pOut) { double *pcf; int rows = n, columns = n + 1; pcf = new double[rows * columns]; if (pcf == nullptr) return 1; //巧妇难为无米之炊，内存都申请不到，还能干嘛！ memcpy(pcf, pCoef, (rows * columns) * sizeof(double)); //据说这个运行效率很高 int r, c, i; //循环变量 int a, b; double x, y; //开始消元，将 pcf 方阵区处理成到直角三角形（直角在右上角）矩阵 for (r = 0; r rows - 1; r++) { //选取主元 a = r; x = _abs(pcf[r * columns + r]); for (i = r + 1; i rows; i++) { //查找主元在哪行 if (x _abs(pcf[i * columns + r])) a = i; } if (a r) { //主元不是当前行（r）,比较麻烦，需要将第 a 行与第 r 行兑换 //第 r 列前面的就不要对换了，因为这些项已经被消元，变成 0 了 for (c = r; c columns; c++) { x = pcf[r * columns + c]; pcf[r * columns + c] = pcf[a * columns + c]; pcf[a * columns + c] = x; } } //开始消元 a = r * columns; //记住将主元的行地址偏移量，以提高程序运行效率 x = -pcf[a + r]; //要多次使用，记下她，以提高程序运行效率 if (x == 0) //主元居然为 0，纯粹是想坑爹，岂能上当！ continue; //继续后面的消元，以便最终判断是无解还是任意解 for (i = r + 1; i rows; i++) { //正在消元 b = i * columns;//记住将要消元的行地址偏移量，以提高程序运行效率 y = pcf[b + r]; //要多次使用，记下她，以提高程序运行效率 if (y != 0) { //y == 0，本行不需要消元 y /= x; //要多次使用，记下她，以提高程序运行效率 pcf[b + r] = 0; //肯定为 0，不用计算。 for (c = r + 1; c columns; c++) pcf[b + c] += pcf[a + c] * y; } } }//至此，pcf 方阵区已经处理成到直角三角形（直角在右上角）矩阵 //回代，将 pcf 方阵区处理成主对角线为 1，其他为 0 的矩阵 int columns_1 = c = columns - 1; //多次用到，提高效率 for (r = rows - 1; r = 1; r--) { b = r * columns; if (pcf[b + r] == 0) { //经过前面的消元，除主元外，其他元应该都为 0 if (pcf[b + columns