【计算方法】三角分解法（线性方程组的求解）

接着我们联立AX=B,A=LU可得到LUX=B;到这里我们假设UX=Y,则原式又转变为LY=B。到这里我们需要逆向思维进行求解，我们采用前向替代法解的Y，然后对UX=Y采用回代法也就是求解上三角线性方程组的方法进行求解，最终解的X。
回代法在之前blog写过了，这里不再赘述。
前向替代法其实就是将回代法反过来用。

A=LU的分解

例题

【问题描述】为求解一个线性方程组，首先采用偏序选主元策略的三角分解法构造矩阵L，U和P，再用前向替换法对方程组LY=PB求解Y，最后用回代法对方程组UX=Y求解X。
【输入形式】在屏幕上依次输入方阵阶数n，系数矩阵A和常数矩阵B。
【输出形式】先输出LU分解结果，再输出方程解。如果系数矩阵A是奇异矩阵，输出: error
【样例1输入】
4
1 2 4 1
2 8 6 4
3 10 8 6
4 12 10 6
21
52
79
82
【样例1输出】
1 0 0 0
0.5 1 0 0
0.25 -0.5 1 0
0.75 0.5 0 1
4 12 10 6
0 2 1 1
0 0 2 0
0 0 0 3
1
2
3
4
【样例1说明】
输入：第1行为方阵阶数4，第2行至5行为系数矩阵A，第6行至9行为常数矩阵B。
输出：第1至第4行输出矩阵L，第5至第8行输出矩阵U，最后4行依次输出方程解：x1, x2, x3, x4。
【评分标准】根据输入得到的输出准确

ACcode：

/*
 * @Author: csc
 * @Date: 2021-04-02 21:12:17
 * @LastEditTime: 2021-04-09 12:45:20
 * @LastEditors: Please set LastEditors
 * @Description: In User Settings Edit
 * @FilePath: \code_formal\course\cal\sanjiaofenjie.cpp
 */
#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>
#define pr printf
#define sc scanf
#define sf(n) scanf("%d", &n)
#define sff(n1, n2) scanf("%d %d", &n1, &n2)
#define sfff(n1, n2, n3) scanf("%d %d %d", &n1, &n2, &n3)
#define sl(n) scanf("%lld", &n)
#define sll(n1, n2) scanf("%lld %lld", &n1, &n2)
#define slll(n1, n2, n3) scanf("%lld %lld %lld", &n1, &n2, &n3)
#define for0(i, n) for (i = 0; i < n; i++)
#define for1n(i, n) for (i = 1; i <= n; i++)
#define foran(i, a, n) for (i = a; i <= n; i++)
#define forna(i, a, n) for (i = n; i >= a; i--)
#define pb push_back
#define fi first
#define se second
#define int long long
#define endl '\n'
#define mem(ara, n) memset(ara, n, sizeof(ara))
#define memb(ara) memset(ara, false, sizeof(ara))
#define all(x) (x).begin(), (x).end()
#define sq(x) ((x) * (x))
#define sz(x) x.size()
const int N = 2e5 + 100;
const int mod = 1e9 + 7;
namespace fastIO
{
   
    inline void input(int &res)
    {
   
        char c = getchar();
        res = 0;
        int f = 1;
        while (!isdigit(c))
        {
   
            f ^= c == '-';
            c = getchar();
        }
        while (isdigit(c))
        {
   
            res = (res << 3) + (res << 1) + (c ^ 48);
            c = getchar();
        }
        res = f ? res : -res;
    }
    inline int qpow(int a, int b)
    {
   
        int ans = 1, base = a;
        while (b)
        {
   
            if (b & 1)
                ans = (ans * base % mod + mod) % mod;
            base = (base * base % mod + mod) % mod;
            b >>= 1;
        }
        return ans;
    }
    int fact(int n)
    {
   
        int res = 1;
        for (int i = 1; i <= n; i++)
        {
   
            res = res * 1ll * i % mod;
        }
        return res;
    }
    int C(int n, int k)
    {
   
        return fact(n) * 1ll * qpow(fact(k), mod - 2) % mod * 1ll * qpow(fact(n - k), mod - 2) % mod;
    }
}
using namespace fastIO;

using namespace std;

int n, i, j;

signed main()
{
   
    int _ = 1;
    //input(_);
    while (_--)
    {
   
        input(n);
        vector<vector<double>> v(n + 1, vector<double>(n + 1));
        for1n(i, n) for1n(j, n) cin >> v[i][j];
        vector<double> b(n + 1);
        for1n(i, n) cin >> b[i];

        // decompose    c 为转置矩阵  v为结果
        vector<vector<double>> c(n + 1, vector<double>(n + 1));
        for1n(i, n) c[i][i] = 1;
        for1n(i, n)
        {
   
            int mai = i;
            foran(j, i, n) if (abs(v[j][i]) > abs(v[mai][i]))
                mai = j;

            if (v[mai][i] == 0) //   eps is enough samll
            {
   
                cout << "error" << endl;
                return 0;
            }

            vector<double> t = v[mai];
            v[mai] = v[i];
            v[i] = t;

            t = c[mai];
            c[mai] = c[i];
            c[i] = t;

            for (j = i + 1; j <= n; j++)
            {
   
                v[j][i] = v[j][i] / v[i][i];  
                for (int k = i + 1; k <= n; k++)
                    v[j][k] = v[j][k] - v[i][k] * v[j][i];
            }
        }
        //  solve
        vector<vector<double>> a = v;

        auto dot = [&](vector<vector<double>> c, vector<double> b) -> vector<double> {
   
            int m = sz(b) - 1;
            vector<double> d(n + 1);
            int k, g;
            for1n(k, m)
                for1n(g, m)
                    d[k] += c[k][g] * b[g];
            return d;
        };
        b = dot(c, b);
        
        vector<double> x(n + 1, 1);
        for1n(i, n)
        {
   
            double m = 0.0;
            foran(j, 1, i - 1)
                m += x[j] * a[i][j];
            x[i] = b[i] - m;
        }

        vector<double> ans(n + 1, 1);
        int k = n;
        while (k > 0)
        {
   
            double m = 0.0;
            foran(j, k + 1, n)
                m = m + ans[j] * a[k][j];
            ans[k] = (x[k] - m) / a[k][k];
            k--;
        }

        // display
        for1n(i, n)
        {
   
            for (j = 1; j <= n; ++j)
                if (j < i)
                    cout << v[i][j] << " ";
                else if (j == i)
                    cout << 1 << " ";
                else
                    cout << 0 << " ";
            cout << endl;
        }
        for1n(i, n)
        {
   
            for (j = 1; j <= n; ++j)
                if (j < i)
                    cout << 0 << " ";
                else
                    cout << v[i][j] << " ";
            cout << endl;
        }

        for1n(i, n) cout << ans[i] << endl;
    }

    return 0;
}