Численные рецепты lubksb рутина не работает

Контекст: я пишу код для решения уравнения Липпмана-Швингера и проверки экспериментальных результатов после использования различных потенциалов.
Чтобы сделать это, я должен вычислить обратную матрицу с разреженной и большой матрицей. Я использую подпрограммы NR ludcmp и lubksb, но что-то не так. После ludcmp матрица имеет много нулей, и lubksb вычисляет вектор b как b [i] / a [i] [j], генерируя -nan большую часть времени. Я не могу понять, что я делаю неправильно.
Спасибо за помощь!

Вот часть моего кода с инверсией матрицы.

    //lu decomp algorithm
void ludcmp(double **a, int n, int *indx, double d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;    //vv stores the implicit scaling of each row.
vv= new double[n];//vector(1,n);
d=1.0; //No row interchanges yet.

for (i=0;i<n;i++)
{//Loop over rows to get the implicit scaling information.
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
//char const *p = "Singular matrix in routine ludcmp"; //
if (big == 0.0) std::cout<<("Singular matrix in routine ludcmp");
//No nonzero largest element.
vv[i]=1.0/big; //Save the scaling.
}

for (j=0;j<n;j++)
{//This is the loop over columns of Crout’s method.
for (i=0;i<j;i++)
{//This is equation (2.3.12) except for i=j.
sum=a[i][j];
for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}

big=0.0;
//Initialize for the search for largest pivot element.
for (i=j;i<n;i++)
{//This is i=j of equation (2.3.12) and i = j+1...N
sum=a[i][j];
for (k=0;k<j-1;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big)
{//Is the figure of merit for the pivot
better than the best so          far?
big=dum;
imax=i;
}
}

if (j != imax)
{//Do we need to interchange rows?
for (k=0;k<n;k++)
{//Yes, do so...
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
d = -d;  //...and change the parity of d.
vv[imax]=vv[j]; //Also interchange the scale factor.
}

indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
/*If the pivot element is zero the matrix is singular
(at least to  the precision of the algorithm).
For some applications on singular matrices,
it is desirable to  substitute TINY for zero.*/
if (j != n)
{//Now, finally, divide by the pivot element.
dum=1.0/(a[j][j]);
for (i=j+1;i<n;i++) a[i][j] *= dum;
}
}
//Go back for the next column in the reduction.

delete[] vv;
}

void lubksb(double **a, int n, int *indx, double b[])
{
int i,ii=0,ip,j;
double sum;
//for (i=1;i<n;i++)
for (i=0;i<n;i++)
{//When ii is set to a positive value,
it will become the index of the first nonvanishing
element of b.
//We now do the forward substitution, equation (2.3.6).
//The only new wrinkle is to unscramble
//the permutation as we go.
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if (ii)
for (j=ii;j<i-1;j++) sum -= a[i][j]*b[j];
else if (sum) ii=i;
//A nonzero element was encountered,
//so from now on we will have to do the sums in the loop above.
b[i]=sum;
}
//for (i=n-2;i>(-1);i--)
for (i=n-1;i>=0;i--)
{//Now we do the backsubstitution, equation (2.3.7).
sum=b[i];
for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
//Store a component of the solution vector X.
}
for(int i=0; i<n; i++) std::cout<<
"b["<<i<<"]= "<<b[i]<<std::endl;
}//matrix inversion

void invmat(double **a, double **y, int N, int *indx, double *b,          double d)
{
ludcmp(a,N,indx,d);
std::ofstream oo; oo.open("invmattttt.txt",std::ios::out);
oo<< "after lu dec" << std::endl;
for(int i=0; i< N; i++)
for(int j=0; j<N; j++) oo << "a["<<i<<"]["<<j<<"]= " << a[i][j] << std::endl;
//Decompose the matrix just once.
//for(int j=1;j<=N;j++)
for(int j=0;j<N;j++)
{
//Find inverse by columns.
//for(int i=1;i<=N;i++) b[i]=0.0;
for(int i=0;i<N;i++) b[i]=0.0;
b[j]=1.0;
lubksb(a,N,indx,b);
oo<< "\n\nafter lubksb" << std::endl;
for(int i=0; i< N; i++)
for(int j=0; j<N; j++) oo <<"a["<<i<<"]["<<j<<"]= " << a[i][j] <<          std::endl;
//for(int i=1;i<N;i++) {y[i][j]=b[i]; std::clog<<i<<std::endl;}          std::clog<<"hey there"<<std::endl;
for(int i=0;i<(N);i++) {y[i][j]=b[i];}}
oo<< "\n\nafter all" << std::endl;
for(int i=0; i< N; i++)
for(int j=0; j<N; j++) oo << y[i][j] << std::endl;
}