运用Doolitle分解法解线性方程组

#include<stdio.h>

const int maxn=1100;

double b[maxn],x[maxn],y[maxn];
double a[1100][1100];
double l[1100][1100];
double u[1100][1100];
int n;

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%lf",&a[i][j]);
    for(int i=1;i<=n;i++) scanf("%lf",&b[i]);
    for(int i=1;i<=n;i++)
        l[i][i]=1;
    for(int k=1;k<=n;k++){
        if(k==1){
            for(int j=1;j<=n;j++) u[1][j]=a[1][j];    
            for(int i=2;i<=n;i++) l[i][1]=a[i][1]/u[1][1];
        }
        else{
            for(int j=k;j<=n;j++){
                double tmp=0;
                for(int s=1;s<=k-1;s++) 
                    tmp+=l[k][s]*u[s][j];
                u[k][j]=a[k][j]-tmp;
            }
            for(int i=k+1;i<=n;i++){
                double tmp=0;
                for(int s=1;s<=k-1;s++)
                    tmp+=l[i][s]*u[s][k];
                l[i][k]=(a[i][k]-tmp)/u[k][k]; 
            }
        }
    }
    y[1]=b[1];
    for(int k=2;k<=n;k++){
        double tmp=0;
        for(int s=1;s<=k-1;s++)
            tmp+=l[k][s]*y[s];
        y[k]=b[k]-tmp;
    }
    x[n]=y[n]/u[n][n];
    for(int k=n-1;k;k--){
        double tmp=0;
        for(int s=k+1;s<=n;s++)
            tmp+=u[k][s]*x[s];
        x[k]=(y[k]-tmp)/u[k][k];
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) printf("%.6f ",l[i][j]);
        puts("");
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) printf("%.6f ",u[i][j]);
        puts("");
    }
    for(int i=1;i<=n;i++) printf("y[%d]=%.6f\n",i,y[i]);
    for(int i=1;i<=n;i++) printf("x[%d]=%.6f\n",i,x[i]);
    return 0;
}
/*
3
1 2 -1
1 -1 5
4 1 -2
3 0 2
*/

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1100;

struct node{
    double m[maxn][maxn];
};

double b[maxn],x[maxn],y[maxn];
node a,l,u;
int n;

int main(){
    cin>>n;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            cin>>a.m[i][j];
    for(int i=1;i<=n;i++) cin>>b[i];
    for(int i=1;i<=n;i++)
        l.m[i][i]=1;
    for(int k=1;k<=n;k++){
        if(k==1){
            for(int j=1;j<=n;j++) u.m[1][j]=a.m[1][j];    
            for(int i=2;i<=n;i++) l.m[i][1]=a.m[i][1]/u.m[1][1];
        }
        else{
            for(int j=k;j<=n;j++){
                double tmp=0;
                for(int s=1;s<=k-1;s++) 
                    tmp+=l.m[k][s]*u.m[s][j];
                u.m[k][j]=a.m[k][j]-tmp;
            }
            for(int i=k+1;i<=n;i++){
                double tmp=0;
                for(int s=1;s<=k-1;s++)
                    tmp+=l.m[i][s]*u.m[s][k];
                l.m[i][k]=(a.m[i][k]-tmp)/u.m[k][k]; 
            }
        }
    }
    y[1]=b[1];
    for(int k=2;k<=n;k++){
        double tmp=0;
        for(int s=1;s<=k-1;s++)
            tmp+=l.m[k][s]*y[s];
        y[k]=b[k]-tmp;
    }
    x[n]=y[n]/u.m[n][n];
    for(int k=n-1;k;k--){
        double tmp=0;
        for(int s=k+1;s<=n;s++)
            tmp+=u.m[k][s]*x[s];
        x[k]=(y[k]-tmp)/u.m[k][k];
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) printf("%.6f ",l.m[i][j]);
        puts("");
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) printf("%.6f ",u.m[i][j]);
        puts("");
    }
    for(int i=1;i<=n;i++) cout<<y[i]<<endl;
    for(int i=1;i<=n;i++) cout<<x[i]<<endl;
    return 0;
}
/*
3
1 2 -1
1 -1 5
4 1 -2
3 0 2
*/