поиск обратной матрицы

[iris_ka]

дана матрица A(почти матрица гильберта-отличается только 1-ая строчка):
1 1 1 1...
1/2 1/3 1/4 1/5...
1/3 1/4 1/5 1/6...
.............................
нужно найти обратную к данной матрице, используя LU- разложение матрицы(разложение на нижнетреугольную и верхнетреуг.)
т.е. в кратце как-то так: для поиска обратной матрицы нужно решить уравнение А*Х=Е(единичная матр.). сначала ищу A=L*U. затем поиск обратной матрицы сводится к решению системы:
L*Y=E
U*X=Y
X-это и есть искомая матрица.

проблема в том, что я не могу найти ошибку в своей программе, очень надеюсь, что найду здесь помощь)
(возможно мой код не блещет изящностью исполнения и экономией памяти, но сейчас главное,чтоб хотя бы работало)
(до строки из "/////////////" все правильно-проверено вручную, на листочке
далее-выкладываю код с комментами, помогите,пожалуйста)

#include <iostream>
#include <math.h>
#include <conio.h>
using namespace std;
 
int main()
{
int i, j, n, k;
cout<<"n=";
cin>>n;   //порядок матрицы
double**matrA=new double*[n];   //массив под матр.А
    for(i=0;i<n;i++)
    {matrA[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {if (i==0) matrA[i][j]=1.0;
    else matrA[i][j]=1.0/(i+j+1);}
    
double**matrL1=new double*[n];   //L1-это вспомогательная нижнетреуг.матр.; при умножении на нее матрицы А,
    for(i=0;i<n;i++)                    //обнуляются поддиагональные элементы А, т.е. L1*A=U
    {matrL1[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {if (i<j) matrL1[i][j]=0.0;
    if (i==j) matrL1[i][j]=1.0;
    if (i>j) matrL1[i][j]=(-1.0)*matrA[i][j]/matrA[j][j];}
    
double**matrU=new double*[n];  //U-верхнетреуг.матр.
    for(i=0;i<n;i++)
    {matrU[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {matrU[i][j]=0.0;
    for (k=0; k<n; k++)
    matrU[i][j]=matrU[i][j]+matrL1[i][k]*matrA[k][j];}
    
double**matrL=new double*[n];  //L-искомая для LU-разложения нижнетреуг.матр.
    for(i=0;i<n;i++)
    {matrL[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {matrL[i][j]=0.0;}
    
double**matrE=new double*[n];  //Е-единичная матр.
    for(i=0;i<n;i++)
    {matrE[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {if (i==j) matrE[i][j]=1.0;
    else matrE[i][j]=0.0;}
//////////////////////////////////    
double sum;                //на данном этапе у нас есть матр.L1 и U,такие, что L1*А=U
for(j=0; j<n; j++)        //хотим "перетащить" L1 в правую часть. ищем нижнетреуг. матрицу L=L1^(-1)
{sum=0.0;                  //т.е. решаем систему L1*L=E прямым ходом
for(i=0; i<n; i++)
{
  for(k=0;k<i;k++) {sum=sum+matrL1[i][k]*matrL[k][j];}
  matrL[i][j]=(matrE[i][j]-sum)/matrL1[i][i];}}
    
double**matrY=new double*[n];   //вспомогательная матр.
    for(i=0;i<n;i++)
    {matrY[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {matrY[i][j]=0.0;}  
      
double**matrX=new double*[n];  // X-искомая матр.(обратная к А)
    for(i=0;i<n;i++)
    {matrX[i]=new double[n];}
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {matrX[i][j]=0.0;}
 
  
for(j=0; j<n; j++)  //решим систему L*Y=E прямым ходом
{sum=0.0;
for(i=0; i<n; i++)
{
  for(k=0;k<i;k++) {sum=sum+matrL[i][k]*matrY[k][j];}
  matrY[i][j]=(matrE[i][j]-sum)/matrL[i][i];}}
  
for(j=n-1; j>=0; j--)    //решим систему U*X=Y обратным ходом
{sum=0.0;
for(i=n-1; i>=0; i--)
{
  for(k=i;k<n;k++) {sum=sum+matrU[i][k]*matrX[k][j];}
  matrX[i][j]=matrY[i][j]-sum;}}    
 
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {cout<<matrX[i][j]<<" ";}
    getchar(); 
    
 return 0;
    
}