Правильно ли решение?(метод Гаусса и LU разложение)

[Сori]

Добрый день!Помогите пожалуйста разобраться, правильно ли решение, которое находит программа(есть сомнения).
Метод Гаусса:

include <iostream>
#include <iomanip>
#include <conio.h>
#include <stdlib.h>
#include <math.h>

int main()
{ /*Введем матрицу размерности 4 на 4 */
	int n=4;
double an[4][4] = { { 2, 7,4,1},
{ 1, 3,5,-2},
{ 1, 5,-9,1},
{5,18,4,5 }};
double bn[4]={ 7,3,1,12}; z =
double xn[4];
	
	for(int k=0;k<n-1;k++)
	{
		for(int i=k+1;i<n;i++)
		{
			double tik;
			tik=an[i][k]/an[k][k];
			bn[i]=bn[i]-an[i][k]*bn[k];
			for(int j=k+1;j<n;j++)
			{
				an[i][j]=an[i][j]-tik*an[k][j];
			};
		};
	};
/*Обратный ход метода Гаусса */

	xn[n-1]=bn[n]/an[n][n];
	double d=0;
	for(k=n-1;k>=0;k--)
	{
		for(int j=k+1;j<n;j++)
		{	d+=an[k][j]*xn[j];
		}
		xn[k]=(bn[k]-d)/an[k][k];
	};
    std::cout<<"X="<<std::endl;
	for(int i=0;i<n;i++)
	{
		std::cout<<std::setw(12)<<xn[i];
		std::cout<<std::endl;

	};
/*Вычисление определителя */
	double det=an[0][0];
	for(int m=1;m<n;m++)
	{
		det=det*an[m][m];
	}
	std::cout<<"detA="<<det<<std::endl;
	getch();
	return 0;
}

Вторая программа(LU разложение):

#include <iostream>
#include <io.h>
#include "stdio.h"
#include "conio.h"
#include <windows.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#define n 4
void main()
{
double T[n][n];
double A[n][n]={{ 15, 7,3,1 },
{ 1, 40, 5,-2 },
{ 1, 5, 15,8 },
{5,18,4,35} 
};
double B[n]={7,3,1,12};
int i;
int j,k;
/*Определим, имеет ли матрица А диагональное преобладание*/
for (i=0;i<n;i++)
{
double p;
p=0;
for(j=0;j<i;j++)
{
P+=A[i][j];
}
for(j=i+1;j<n;j++)
{
P+=A[i][j];
}
if(fabs(A[i][i])<p)
{
std::cout<<"net diagonalnogo preobladaniya"<<std::endl;
getch();
}
else
{	
std::cout<<"diagonalnoe preobladanie"<<std::endl;
}
/*Создаем матрицу LU-разложения – T.*/
for (i=1;i<=n;i++)
{
T[1][i]=A[1][i]; 
for(i=2;i<=n;i++)
{
T[i][1]=A[i][1]/T[1][1]; 
}
}
for(i=2;i<=n;i++)
{
for(j=i;j<=n;j++)
{
p=A[i][j];
for(k=1;k<i;k++)
{
p-=T[i][k]*T[k][j];
}
T[i][j]=p;
}
for(j=i+1;j<=n;j++)
{
p=A[j][i];}
for(k=1;k<i;k++);
{
p-=T[j][k]*T[k][i];
}
T[j][i]=p/T[i][i];
}
std::cout<<"T="<<std::endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
/*Вывод матрицы Т*/
std::cout<<T[i][j];
}
std::cout<<std::endl;
};

for(i=1;i<n;i++)
{
p=B[i];
for(j=0;j<i;j++)
{
p-=T[i][j]*B[j];
}
B[i]=p;
}
std::cout<<"Y="<<std::endl;
for(i=0;i<n;i++)
{
std::cout<<B[i];
std::cout<<std::endl;
};
double x[4];
x[n-1]=B[n-1]/T[n-1][n-1];
double t=0;
for(k=n-2;k>=0;k--)
{
for(j=k+1;j<n;j++)
{	
t+=T[k][j]*x[j];
}
x[k]=(B[k]-t)/T[k][k];
};
std::cout<<"X="<<std::endl;
for(i=0;i<n;i++)
{
std::cout<<x[i];
std::cout<<std::endl;
};
/*Вычисление определителя матрицы А*/
double det=T[0][0];
int f
for(f=1;f<n;f++)
{
det=det*T[f][f];
}
std::cout<<"detA="<<det<<std::endl;
getch();
}
}