高斯消元法

高斯消元法

（以下摘自百度百科：高斯消元法-百度百科）
内容
消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代人到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。
核心
1)两方程互换，解不变；
2)一方程乘以非零数k，解不变；
3)一方程乘以数k加上另一方程，解不变

找出逆矩阵

高斯消元法可以用来找出一个可逆矩阵的逆矩阵。设A 为一个N * N的矩阵，其逆矩阵可被两个分块矩阵表示出来。将一个N * N单位矩阵 放在A 的右手边，形成一个N * 2N的分块矩阵B = [A,I] 。经过高斯消元法的计算程序后，矩阵B 的左手边会变成一个单位矩阵I ，而逆矩阵A ^(-1) 会出现在B 的右手边。
假如高斯消元法不能将A 化为三角形的格式，那就代表A 是一个不可逆的矩阵。
应用上，高斯消元法极少被用来求出逆矩阵。高斯消元法通常只为线性方程组求解。

分析

高斯消元法的算法复杂度是O(n^3)；这就是说，如果系数矩阵的是n × n，那么高斯消元法所需要的计算量大约与
成比例。
高斯消元法可用在任何域中。
高斯消元法对于一些矩阵来说是稳定的。对于普遍的矩阵来说，高斯消元法在应用上通常也是稳定的，不过亦有例外。

推荐一篇比较好理解的blog ：高斯消元法
代码：
（可增加一列变成增广矩阵。。。。然后回代求解。。。。）

//#include<bits/stdc++.h>
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<stack>
#include<queue>			//stack的头文件
using namespace std ;
typedef long long ll;
#define MAXN 105
#define INF 0x3f3f3f3f
#define MALL (BiTnode *)malloc(sizeof(BiTnode));

double a[MAXN][MAXN];
void init(int n)						//输入n x n的矩阵；
{
    memset(a, 0, sizeof(a));
    for(int i=1; i<=n; ++i)
        for(int j=1; j<=n; ++j)
            scanf("%lf", &a[i][j]);
}

int gaosi(int n)					//高斯消元（上三角）
{
    for(int i=1; i<n; ++i)
    {
        if(!a[i][i])
            return 0;
        for(int j=i+1; j<=n; ++j)
        {
            double ans = 1.0*a[j][i]/a[i][i];
            for(int k=i; k<=n; ++k)
                a[j][k] = a[j][k] - 1.0*ans*a[i][k];
        }
    }
    return 1;
}

void print(int n)
{
    for(int i=1; i<=n; ++i)
    {
        for(int j=1; j<=n; ++j)
            printf("%.2f ", a[i][j]);
        cout << '\n';
    }
}

int main()
{
    int n;
    cin >> n;
    init(n);
    int i = gaosi(n);
    print(n);
    return 0;
}

最近翻到了一篇很好的模板，特此放上来：高斯消元快速入门

这一篇的模板分整型，double，异或，活脱脱的大牛：高斯消元模板

 #include<stdio.h>
    #include<algorithm>
    #include<iostream>
    #include<string.h>
    #include<math.h>
    using namespace std;
    const int MAXN=50;
    int a[MAXN][MAXN];//增广矩阵
    int x[MAXN];//解集
    bool free_x[MAXN];//标记是否是不确定的变元
    int gcd(int a,int b){
        if(b == 0) return a; else return gcd(b,a%b);
    }
    inline int lcm(int a,int b){
        return a/gcd(a,b)*b;//先除后乘防溢出
    }
    // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，
    //-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
    //有equ个方程，var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
    int Gauss(int equ,int var){
        int i,j,k;
        int max_r;// 当前这列绝对值最大的行.
        int col;//当前处理的列
        int ta,tb;
        int LCM;
        int temp;
        int free_x_num;
        int free_index;

        for(int i=0;i<=var;i++){
            x[i]=0;
            free_x[i]=true;
        }

        //转换为阶梯阵.
        col=0; // 当前处理的列
        for(k = 0;k < equ && col < var;k++,col++){// 枚举当前处理的行.
        // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
            max_r=k;
            for(i=k+1;i<equ;i++){
                if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
            }
            if(max_r!=k){// 与第k行交换.
                for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
            }
            if(a[k][col]==0){// 说明该col列第k行以下全是0了，则处理当前行的下一列.
                k--;
                continue;
            }
            for(i=k+1;i<equ;i++){// 枚举要删去的行.
                if(a[i][col]!=0){
                    LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                    ta = LCM/abs(a[i][col]);
                    tb = LCM/abs(a[k][col]);
                    if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
                    for(j=col;j<var+1;j++){
                        a[i][j] = a[i][j]*ta-a[k][j]*tb;
                    }
                }
            }
        }
        // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
        for (i = k; i < equ; i++){ // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
            if (a[i][col] != 0) return -1;
        }
        // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
        // 且出现的行数即为自由变元的个数.
        if (k < var){
            return var - k; // 自由变元有var - k个.
        }
        // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
        // 计算出Xn-1, Xn-2 ... X0.
        for (i = var - 1; i >= 0; i--){
            temp = a[i][var];
            for (j = i + 1; j < var; j++){
                if (a[i][j] != 0) temp -= a[i][j] * x[j];
            }
            if (temp % a[i][i] != 0) return -2; // 说明有浮点数解，但无整数解.
            x[i] = temp / a[i][i];
        }
        return 0;
    }
    int main(void){
    //    freopen("in.txt", "r", stdin);
    //    freopen("out.txt","w",stdout);
        int i, j;
        int equ,var;
        while (scanf("%d %d", &equ, &var) != EOF){
            memset(a, 0, sizeof(a));
            for (i = 0; i < equ; i++){
                for (j = 0; j < var + 1; j++){
                    scanf("%d", &a[i][j]);
                }
            }
            int free_num = Gauss(equ,var);
            if (free_num == -1) printf("无解!\n");
            else if (free_num == -2) printf("有浮点数解，无整数解!\n");
            else if (free_num > 0){
                printf("无穷多解! 自由变元个数为%d\n", free_num);
                for (i = 0; i < var; i++){
                    if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
                    else printf("x%d: %d\n", i + 1, x[i]);
                }
            }else{
                for (i = 0; i < var; i++){
                    printf("x%d: %d\n", i + 1, x[i]);
                }
            }
            printf("\n");
        }
        return 0;
    }