数值实验-高斯消元LU分解

实验课的作业，用LU分解矩阵A求解线性方程组。

注：只是课程的设计，不可以当作ACM的模板，也是不是考虑的很周到，只是提供一个思路。

#include "cstdio"
#include "cstring"
#include "cstdlib"
#include "cmath"
using namespace std;

const double eps = 1e-8;
const int maxr = 100;
const int maxc = 100;

struct Matrix {
    double m[maxr][maxc];
    int row, col;
};

Matrix operator - (Matrix a, Matrix b) {
    Matrix c;
    c.col = a.col; c.row = a.row;
    for (int i = 0; i < a.row; i++) {
        for (int j = 0; j < a.col; j++) {
            c.m[i][j] = a.m[i][j]-b.m[i][j];
        }
    }
    return c;
}

Matrix operator + (Matrix a, Matrix b) {
    Matrix c;
    c.col = a.col; c.row = a.row;
    for (int i = 0; i < a.row; i++) {
        for (int j = 0; j < a.col; j++) {
            c.m[i][j] = a.m[i][j] + b.m[i][j];
        }
    }
    return c;
}

Matrix operator * (Matrix a, Matrix b) {
    Matrix c;
    memset(c.m, 0, sizeof(c.m));
    c.row = a.row; c.col = b.col;
    for (int i = 0; i < a.row; i++) {
        for (int k = 0; k < a.col; k++) {
            for (int j = 0; j < b.col; j++) {
                c.m[i][j] += a.m[i][k]*b.m[k][j];
            }
        }
    }
    return c;
}

Matrix pre_solve(Matrix a, Matrix b) {
    Matrix c;
    c.row = b.row; c.col = 1;
    for (int i = 0; i < a.row; i++) {
        double t = b.m[i][0];
        for (int j = 0; j < i; j++) {
            t -= a.m[i][j]*c.m[j][0];
        }
        c.m[i][0] = t/a.m[i][i];
    }
    return c;
}
//后代法
Matrix pos_solve(Matrix a, Matrix b) {
    Matrix c;
    c.row = b.row; c.col = 1;
    for (int i = a.row-1; i >= 0; i--) {
        double t = b.m[i][0];
        for (int j = i+1; j < a.col; j++) {
            t -= a.m[i][j]*c.m[j][0];
        }
        c.m[i][0] = t/a.m[i][i];
    }
    return c;
}
//生成初等行变换矩阵
Matrix init_Matrix(int p, int q, int n) {
    Matrix c;
    c.col = c.row = n;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (i == j) c.m[i][j] = 1.0;
            else c.m[i][j] = 0.0;
        }
    }
    if (p == q) return c;
    c.m[p][p] = 0.0;
    c.m[p][q] = 1.0;
    c.m[q][q] = 0.0;
    c.m[q][p] = 1.0;
    return c;
}

void show_Matrix(Matrix a) {
    for (int i = 0; i < a.row; i++) {
        printf("%.8f", a.m[i][0]);
        for (int j = 1; j < a.col; j++) {
            printf(" %.8f", a.m[i][j]);
        }
        printf("\n");
    }
}
//高斯消元
Matrix Gauss(Matrix a, Matrix b) {
    Matrix c, L, U = a;
    for (int i = 0; i < a.col-1; i++) {
        for (int j = i+1; j < a.row; j++) {
            U.m[j][i] = U.m[j][i]/U.m[i][i];
        }
        for (int j = i+1; j < a.row; j++) {
            for (int k = i+1; k < a.col; k++) {
                U.m[j][k] -= U.m[i][k]*U.m[j][i];
            }
        }
    }
    L.col = U.col;L.row = U.row;
    for (int i = 0; i < U.row; i++) {
        for (int j = 0; j < U.col; j++) {
            if (i == j) L.m[i][j] = 1.0;
            else if (i > j) L.m[i][j] = U.m[i][j], U.m[i][j] = 0.0;
            else L.m[i][j] = 0.0;
        }
    }
    c = pre_solve(L, b);
    return pos_solve(U, c);
}
//高斯列主元
Matrix ad_Gauss(Matrix a, Matrix b) {
    Matrix c, L, U = a, P = init_Matrix(0, 0, a.row);
    L.col = U.col; L.row = U.row;
    memset(L.m, 0, sizeof(L.m));
    for (int i = 0; i < a.col-1; i++) {
        int K = i;
        double maxx = 0.00;
        for (int j = i; j < a.row; j++) {
            if (fabs(U.m[i][j]) - maxx > eps) {
                maxx = fabs(U.m[i][j]), K = j;
            }
        }
        P = init_Matrix(i, K, a.row)*P;
        U = init_Matrix(i, K, a.row)*U;
        for (int j = i+1; j < U.row; j++) {
            U.m[j][i] = U.m[j][i]/U.m[i][i];
        }
        for (int j = i+1; j < U.row; j++) {
            for (int k = i+1; k < U.col; k++) {
                U.m[j][k] -= U.m[i][k]*U.m[j][i];
            }
        }
    }
    L.col = U.col;L.row = U.row;
    for (int i = 0; i < U.row; i++) {
        for (int j = 0; j < U.col; j++) {
            if (i == j) L.m[i][j] = 1.0;
            else if (i > j) L.m[i][j] = U.m[i][j], U.m[i][j] = 0.0;
            else L.m[i][j] = 0.0;
        }
    }
    c = pre_solve(L, P*b);
    return pos_solve(U, c);
}

int main(int argc, char const *argv[])
{
    Matrix A, B, C;
    int n;
    printf("输入系数矩阵的维数\n");
    scanf("%d", &n);
    A.col = A.row = n;
    B.row = n; B.col = 1;
    printf("输入矩阵A\n");
    for (int i = 0; i < A.row; i++) {
        for (int j = 0; j < A.col; j++) scanf("%lf", &A.m[i][j]);
    }
    printf("输入矩阵B\n");
    for (int i = 0; i < B.row; i++) scanf("%lf", &B.m[i][0]);
    C = ad_Gauss(A, B);
    printf("利用高斯消元得到的解为\n");
    show_Matrix(C);
    return 0;
}