数值分析2.2 直接三角分解方法

#include <iostream>
#include <cstdio>
#include <math.h>
using namespace std;
const int MAXN = 1000;
double a[MAXN][MAXN];
double l[MAXN][MAXN];
double u[MAXN][MAXN];
double b[MAXN];
double x[MAXN];
double y[MAXN];
int main()
{
    int n;
    double maxs;
    freopen("in.txt","r",stdin);
    while(cin >> n){
        for(int i = 0;i<n;i++)cin >> b[i];
        for(int i = 0;i<n;i++){
            for(int j = 0;j<n;j++){
                cin >> a[i][j];
            }
        }
        //memset(l,0,sizeof(l));
        //memset(u,0,sizeof(u));
        for(int i = 0;i<n;i++)l[i][i] = 1;
        for(int k = 0;k<n;k++){///求第k行u和第k列l
            for(int j = k;j<n;j++){
                    u[k][j] = a[k][j];
                for(int m = 0;m<k;m++){
                    u[k][j]-=l[k][m]*u[m][j];
                }
            }
            for(int j = k+1;j<n;j++){
                    l[j][k] = a[j][k];
                for(int m = 0;m<k;m++){
                    l[j][k]-=l[j][m]*u[m][k];
                }
                l[j][k]/=u[k][k];
            }
        }
        ///回带求y
        for(int i = 0;i<n;i++){
                y[i] = b[i];
            for(int j = 0;j<i;j++){
                y[i]-=l[i][j]*y[j];
            }
        }
        ///回带求x
        for(int i = n-1;i>=0;i--){
                x[i] = y[i];
            for(int j = n-1;j>i;j--){
                x[i]-=u[i][j]*x[j];
            }
            x[i]/=u[i][i];
        }
        ///输出
        for(int i = 0;i<n;i++){
            if(!i)cout << x[i];
            else cout << " " << x[i];
        }
        cout << endl;
    }
    return 0;
}