用高斯消元法解三元一次方程组，C语言

参阅我的文章：http://wenku.baidu.com/view/d4ea2273650e52ea5418981d.html

#include "stdafx.h" //VS 预编译头文件，其他系统请删除
#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<math.h>
#include<time.h>

//VS 2013 否决了 scanf 等函数，为了使用，加上下句。
//其他系统请删除
#pragma warning(disable:4996)
int GaussJordanElimination(int n, const double *pCoef, double *pOut);
//VS 主函数签名格式。其他系统请改变签名，如：
//int main()
int _tmain(int argc, _TCHAR* argv[])
{
double cf[3][4] = { {-0.02, 2.0, 2.0, 0.4}, {1.0, 0.78125, 0.0, 1.3816}, {3.996, 5.526, 4.0, 7.4178} };
double rs[3];
int i;
i = GaussJordanElimination(3, (double*)cf, rs);
printf("x1 = %lf, x2 = %lf, x3 = %lf\n", rs[0], rs[1], rs[2]);
system("pause"); //避免窗口一闪而退
return 0;
}
//绝对值函数
__inline double _abs(double v)
{
return v < 0 ? -v : v;
}

//线性方程组列主元高斯消元法
//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];
pcf = (double*)malloc(rows * columns * sizeof(double));
if (pcf == 0) 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_1] == 0)
{ //常数项为 0，方程有不定解
free(pcf);
return 2;
}
else
{ //常数项为 0，方程有无解
free(pcf); //释放内存
return -1;
}
}
pcf[b + columns_1] /= pcf[b + r];
pcf[b + r] = 1; //肯定为 1，不用计算。
y = -pcf[b + columns_1];
//回代
for (i = r - 1; i >= 0; i--)
{
pcf[i * columns + columns_1] += pcf[i * columns + r] * y;
pcf[i * columns + r] = 0; //已经回代，此项已消，置为 0。
}
}
//处理第一行数据
pcf[columns_1] /= pcf[0];
pcf[0] = 1;
//至此，回代过程结束，pcf 矩阵的最后一列就是结果
//返回结果
for (r = 0; r < rows; r++)
pOut[r] = pcf[r * columns + columns_1];
free(pcf);
return 0;
}本回答被提问者和网友采纳