高斯消元法【模板】

高斯消元法，消成行阶梯型矩阵。

下面两种消元法的时间复杂度都是 $O(n^3)$.

#include<bits/stdc++.h>
using namespace std;

const int maxn = 100+10;
typedef double Matrix[maxn][maxn];

//要求系数矩阵可逆
//这里的A是增广矩阵，即A[i][n] 是第i个方程右边的常数bi
//运行结束后A[i][n] 是第i个未知数的值
void gauss_elimination(Matrix A, int n)
{
    int i, j, k, r;

    for(i = 0;i < n;i++)  //消元过程
    {
        //选绝对值一行r并与第i行交换
        r = i;
        for(j = i+1; j < n;j++)
            if(fabs(A[j][i] > fabs(A[r][i]))) r = j;
        if(r != i) for(j = 0;j <= n;j++)  swap(A[r][j], A[i][j]);

        //与第i+1~n行进行消元
        for(k = i+1; k < n;k++)
        {
            double f = A[k][i] / A[i][i];
            for(int j = i;j <= n;j++)  A[k][j] -= f * A[i][j];      //已经是阶梯型矩阵了，所以从i开始
        }
    }

    //回代过程
    for(i = n-1;i >= 0;i--)
    {
        for(j = i+1; j < n;j++)
            A[i][n] -= A[j][n] * A[i][j];
        A[i][n] /= A[i][i];
    }
}

int n;
Matrix M;

int main()
{
    while(scanf("%d", &n) == 1)
    {
        for(int i = 0;i < n;i++)
            for(int j = 0;j <= n;j++)
                scanf("%lf", &M[i][j]);
        gauss_elimination(M, n);

        for(int i = 0;i < n;i++)  printf("%.8f
", M[i][n]);
    }
}

高斯-约当消元法，消成对角矩阵，从而省略掉回代过程。

#include<bits/stdc++.h>
using namespace std;

const double eps = 1e-8;
const int maxn = 100+10;
typedef double Matrix[maxn][maxn];

//结果为A[i][n]/A[i][i]
void gauss_jordan(Matrix A, int n)
{
    int i, j, k, r;
    for(i = 0;i < n;i++)
    {
         //选绝对值一行r并与第i行交换
        r = i;
        for(j = i+1;j < n;j++)
            if(fabs(A[j][i]) > fabs(A[r][i]))  r = j;
        if(fabs(A[r][i]) < eps)  continue;      //放弃这一行，直接处理下一行
        if(r != i)  for(j = 0;j <= n;j++)  swap(A[r][j], A[i][j]);

        //与除第i行外的其他行进行消元
        for(k = 0;k < n;k++)  if(k != i)
            for(j = n;j >= i;j--)  A[k][j] -= A[k][i] / A[i][i] * A[i][j];
    }
}


int n;
Matrix M;

int main()
{
    while(scanf("%d", &n) == 1)
    {
        for(int i = 0;i < n;i++)
            for(int j = 0;j <= n;j++)
                scanf("%lf", &M[i][j]);
        gauss_jordan(M, n);

        for(int i = 0;i < n;i++)  printf("%.8f
", M[i][n] / M[i][i]);
    }
}

Code From:

《算法竞赛入门经典训练指南》——刘汝佳、陈锋编著