LU-разложение

#include "stdafx.h"
#include "stdio.h"
#include "conio.h"
#include "math.h"
#include "iostream"
//using namespace std;
 
int const N=5;
double a[N][N],l[N][N],u[N][N];
double b[N],y[N],x[N];
int s,i,j,n;
double sum,tmp;
 
int main(){
 for(i=0;i<N;i++){
     for(j=0;j<N;j++){
         l[i][j]=0;
         u[i][j]=0;
     }
     y[i]=0;
     x[i]=0;
     u[i][i]=1;
 }
 //==========================================================
 cin>>n;
 for(i=0;i<n;i++)
  for(j=0;j<n;j++)
   cin>>a[i][j]; 
 for(i=0;i<n;i++) 
   cin>>b[i];
//===========================================================
 for(i=0;i<n;i++)
     for(j=0;j<n;j++){
         sum=0;
         if(i>=j){
             for(s=0;s<=j-1;s++)
                 sum+=(l[i][s])*(u[s][j]);
             l[i][j]=a[i][j]-sum;
         } else {
             for(s=0;s<=i-1;s++)
                 sum+=l[i][s]*u[s][j];
             u[i][j]=(a[i][j]-sum)/l[i][i];
         }
     }
//======================================================================================
 y[0]=b[0]/l[0][0];
 for(i=1;i<n;i++){
     sum=0;
     for(j=0;j<i;j++)
         sum+=l[i][j]*y[j];
     y[i]=(b[i]-sum)/l[i][i];
 }
 
 x[n-1]=y[n-1];
 for(i=n-2;i>=0;i--){
     sum=0;
     for(j=n-1;j>=i+1;j--)
         sum+=x[j]*u[i][j];
     x[i]=y[i]-sum;
 }
 
 for(i=0;i<n;i++)
     cout<<x[i]<<' ';
 cout<<'\n';
 for(i=0;i<n;i++){
     for(j=0;j<n;j++)
         cout<<l[i][j]<<' ';
     cout<<endl;
     }
 for(i=0;i<n;i++){
     for(j=0;j<n;j++)
         cout<<u[i][j]<<' ';
     cout<<endl;
     }
 
     printf("%4lf ",x[i]);
 getch();
 return 0;
}