计算方法 - 矩阵的三角分解 - 解三对角线方程组的追赶法（线性方程组的解法）

详细过程见《数值计算方法》
题目：

书上的LU矩阵：

最开始用书上写的，写完之后发现老师样例要求的跟书上是反着的，所以无奈只能重新推导，再打了遍代码，实属是无语

题目要求版本及推导：

书上版本C++代码：

#include <iostream>
#include <cstdio>
#include <queue>
#include <cstring>
#include <cmath>
using namespace std;
#define ll long long
const int maxn = 205;
const int INF = 0x3f3f3f;
double temp[105 * 105] = {0.0};
double x[105] = {0.0}, y[105] = {0.0};
double Alpha[1005] = {0.0}, Beta[1005] = {0.0}, Gamma[1005] = {0.0};
int n = 0, n2 = 0;

inline void IN() { //输入
    int t = 0;
    while(scanf("%lf",&temp[++t]) != EOF) {}
    while(n*(n+1) != (t-1)) {n++;}
    n2 = n * n;
}

inline double *A(int i, int j) {
    return &temp[j + (i-1) * n];
}

inline double *f(int j) {
    return &temp[j + n2];
}

inline double *b(int i) {
    return &temp[i + (i-1) * n];
}

inline double *c(int i) {
    return &temp[i + 1 + (i-1) * n];
}

inline double *a(int i) {
    return &temp[i - 1 + (i-1) * n];
}

inline void OUT() {
    printf("L:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(i == j) printf("%.5lf ", Alpha[i]);
            else if(j == i-1) printf("%.5lf ", Gamma[i]);
            else printf("0.00000 ");
        }
        printf("\n");
    }
    printf("U:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(i == j) printf("1.00000 ");
            else if(j == i+1) printf("%.5lf ", Beta[i]);
            else printf("0.00000 ");
        }
        printf("\n");
    }
    printf("y:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", y[i]);
    }
    printf("x:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", x[i]);
    }
}

inline void get_alpha_beta_gamma(){
    Alpha[1] = *b(1);
    for(int i = 2; i <= n; ++i) {
        Gamma[i] = *a(i);
        Beta[i-1] = *c(i-1) / Alpha[i-1];
        Alpha[i] = *b(i) - *a(i) * Beta[i-1]; 
    }
}

inline void get_xy() {
    y[1] = *f(1) / *b(1);
    for(int i = 1; i <= n; ++i) {
        y[i] = (*f(i) - *a(i) * y[i-1]) / (*b(i) - *a(i) * Beta[i-1]);
    }
    x[n] = y[n];
    for(int i = n-1; i >= 1; --i) {
        x[i] = y[i] - Beta[i] * x[i+1];
    }
}

int main() {
    IN();
    get_alpha_beta_gamma();
    get_xy();
    OUT();
    return 0;
}

输入输出结果：

题目版本C++代码：

#include <iostream>
#include <cstdio>
#include <queue>
#include <cstring>
#include <cmath>
using namespace std;
#define ll long long
const int maxn = 205;
const int INF = 0x3f3f3f;
double temp[105 * 105] = {0.0};
double x[105] = {0.0}, y[105] = {0.0};
double Alpha[1005] = {0.0}, Beta[1005] = {0.0}, Gamma[1005] = {0.0};
int n = 0, n2 = 0;

inline void IN() { //输入
    int t = 0;
    while(scanf("%lf",&temp[++t]) != EOF) {}
    while(n*(n+1) != (t-1)) {n++;}
    n2 = n * n;
}

inline double *A(int i, int j) {
    return &temp[j + (i-1) * n];
}

inline double *f(int j) {
    return &temp[j + n2];
}

inline double *b(int i) {
    return &temp[i + (i-1) * n];
}

inline double *c(int i) {
    return &temp[i + 1 + (i-1) * n];
}

inline double *a(int i) {
    return &temp[i - 1 + (i-1) * n];
}

inline void OUT() {
    printf("L:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(i == j) printf("1.00000 ");
            else if(j == i-1) printf("%.5lf ", Gamma[i]);
            else printf("0.00000 ");
        }
        printf("\n");
    }
    printf("U:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(i == j) printf("%.5lf ", Alpha[i]);
            else if(j == i+1) printf("%.5lf ", Beta[i]);
            else printf("0.00000 ");
        }
        printf("\n");
    }
    printf("y:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", y[i]);
    }
    printf("x:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", x[i]);
    }
}

inline void get_alpha_beta_gamma(){
    Alpha[1] = *b(1); Beta[1] = *c(1);
    for(int i = 2; i <= n; ++i) {
        Beta[i] = *c(i);
        Gamma[i] = *a(i) / Alpha[i-1];
        Alpha[i] = *b(i) - Gamma[i] * Beta[i-1]; 
    }
}

inline void get_xy() {
    y[1] = *f(1);
    for(int i = 2; i <= n; ++i) {
        y[i] = *f(i) - Gamma[i] * y[i-1];
    }
    x[n] = y[n] / Alpha[n];
    for(int i = n-1; i >= 1; --i) {
        x[i] = (y[i] - Beta[i] * x[i+1]) / Alpha[i];
    }
}

int main() {
    IN();
    get_alpha_beta_gamma();
    get_xy();
    OUT();
    return 0;
}

输入输出结果：