LU分解でC言語

#連立方程式をLU分解で解いたのでC言語のソースだけ貼っておく

#include<iostream>
#include<stdio.h>
#define N 3
/*
2 5 7 4
4 13 20 11
8 29 50 29
*/
using namespace std;
int main(){
    double mat[N][N+1]={{2,5,7,4},{4,13,20,11},{8,29,50,29}};//元々の拡大係数行列
    double mat_l[N][N]={{0,},{0,}}; //L
    double mat_u[N][N]={{0,},{0,}};//U
    double pivot = 0.0,tmp=0.0;    
    double fl[N]={4,11,29};//拡大係数成分格納 
    double u[N]={0,};//行列uの解
    double l[N]={0,};//最終解
    int i=0,j=0;//カウント変数
    for(int i=0;i<N;i++)//Lを単位行列にする
        mat_l[i][i]=1;

    
    for( i = 0; i < N; i++){
        for( j = 0; j < N; j++){
            mat_u[i][j]=mat[i][j];//行列コピー
        }
    }
    
    for(i = 0; i < N ; i++){
        pivot=mat_u[i][i];//pivotを取り出す
        for(j=i+1;j<N;j++){
            tmp=mat_u[j][i];//被掃き出し行のpivot列の成分を取り出してtmpに代入
            for(int k=0; k < N ;k++){
                mat_u[j][k] = mat_u[j][k] - (tmp/pivot)*mat_u[i][k];//行列Uを掃き出し法で上三角行列にする
             }
                mat_l[j][i]=tmp/pivot;//行列Lの被掃き出し行,pivot列の成分にtmpをpivotで割ったものを代入すると
                                     //なぜかちゃんとLU分解の行列Lになる
        }
    }    
    u[0]=mat_l[0][0]*fl[0];

    for(i=1;i<N;i++){
        u[i]=fl[i];//拡大係数成分を代入
        for(j=0;j<i;j++){ 
         u[i]=u[i]-mat_l[i][j]*u[j];//行列Lの解を求める(代入法）
        }
    
    }

    l[N-1]=mat_u[N-1][N-1]/u[N-1];//zを求める（N*N+1)拡大係数行列なら最終行の変数の解

    for(i=N-2;i >= 0 ;i--){
        l[i]=u[i];
        for(j=N-1; j > i;j-- ){
            
           l[i]=l[i] - mat_u[i][j]*l[j];
        }
        l[i]=l[i]/mat_u[i][i];
    }

  

  for( i = 0; i < 3; i++)
  {
      cout << l[i] << endl;
  }
    
return 0;
}

クッソ適当です