Матрица обратная между Matlab / Octave и переполнением стека

РЕДАКТИРОВАТЬ: извините плохой пример, я пытался получить те же цифры, что и в моей программе, но не повезло, я не могу получить пример для работы.
Единственное, что я могу сделать, это показать вам разницу между Octave / Matlab и C ++

например, в Octave / Matlab у меня есть: 3.8070e + 010 и в C ++: 38121149284.106712

или -1,4971e + 011 в ML / O и -149962686456,46307 в C ++

так что вы видите, что есть разница.

Пожалуйста, дайте мне несколько советов, где искать другой подход.

функция А:

void InvdiagMat(double **Mat, int NbElement)
{
double ** temp;
int i;
temp=new double*[NbElement];
for(i=0;i<NbElement;i++)
{
temp[i]=new double [NbElement];
for(int j=0;j<NbElement;j++)
{
temp[i][j]=0;
}
}

for(i=0;i<NbElement;i++)
{

for(int j=0;j<NbElement;j++)
{
temp[i][i]=1/Mat[i][i];
if(j!=i)
{
temp[i][j]=-Mat[i][j]/Mat[i][i];
}
for(int k=0;k<NbElement;k++)
{

if(k!=i)
{
temp[k][i]=Mat[k][i]/Mat[i][i];
}
if(j!=i &&k!=i)
{
temp[k][j]=Mat[k][j]-Mat[i][j]*Mat[k][i]/Mat[i][i];
}
}

}
for(int i=0;i<NbElement;i++)
{
for(int j=0;j<NbElement;j++)
{
Mat[i][j]=temp[i][j];
}

}
}

}

а также функция B:

    double** invdet(double** a, int n)
/* This function computes both the determinant of matrix a and its inverse matrix */

{
int i,j,k,l,m,*indx;
double d,*col;

double ** inv;
inv = new double*[n];
for (j=0; j<n; j++)
{
inv[j] = new double[n];
}

col=vector(0, n-1);
indx=ivector(0, n-1);

ludcmp(a, n, indx, &d);
for (j=0; j<n; j++)
{
d *= a[j][j];
for (i=0; i<n; i++)
col[i]=0.0;
col[j]=1.0;
lubksb(a,n,indx,col);
for (i=0; i<n; i++)
inv[i][j]=col[i];
}

return inv;
}void ludcmp(double** aa, int n, int *indx, double*d)
{
int i,imax=0,j,k;
double big,dum,sum,temp;
double *vv;

vv=vector(0,n-1);
*d=1.0;
for (i=0; i<n;i++) {
big=0.0;
for (j=0; j<n; j++)
if ((temp=fabs(aa[i][j])) > big) big=temp;
//if (big == 0.0) printf("Singular matrix\n");
vv[i]=1.0/big;
}
for (j=0; j<n; j++)
{
for (i=0; i<j; i++) {
sum=aa[i][j];
for (k=0; k<i; k++)
sum -= aa[i][k]*aa[k][j];
aa[i][j]=sum;
}
big=0.0;
for (i=j; i<n; i++)
{
sum=aa[i][j];
for (k=0; k<j; k++)
sum -= aa[i][k]*aa[k][j];
aa[i][j]=sum;
if ((dum = vv[i] * fabs(sum)) >= big)
{
big = dum;
imax = i;
}
}
if (j != imax)
{
for (k=0; k<n; k++)
{
dum = aa[imax][k];
aa[imax][k] = aa[j][k];
aa[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j]=imax;
if (aa[j][j] == 0.0)
aa[j][j]=TINY;
if (j != n-1)
{
dum=1.0/aa[j][j];
for (i = j+1; i < n; i++)
aa[i][j] *= dum;
}
}
}

void lubksb(double **a,int n,int *indx,double b[])
{
int i,ii=0,ip,j;
double sum;

for (i=0; i<n; i++) {
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if (ii>=0)
for (j=ii; j<=i-1; j++) sum -= a[i][j]*b[j];
else if (sum) ii=i;
b[i]=sum;
}
for (i=n-1; i>=0; i--) {
sum=b[i];
for (j=i+1; j<n; j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
}
}int *ivector(int nl,int nh)
{
int *v;
v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
if (!v) printf("allocation failure in ivector()");
return v-nl;
}double *vector(int nl,int nh)
{
double *v;
v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
if (!v) printf("allocation failure in vector()");
return v-nl;
}void free_vector(double *v,int nl,int nh)
{
free((char*) (v+nl));
}

double **matrix(int nrl,int nrh,int ncl,int nch)
{
int i;
double **m;
m=(double **)malloc((unsigned) (nch-ncl+1)*sizeof(double*));
if (!m) printf("allocation failure 1 in matrix()");
m -= nrl;

for (i=nrl; i<=nrh; i++) {
m[i]=(double *)malloc((unsigned) (nch-ncl+1)*sizeof(double));
if (!m[i]) printf("allocation failure 2 in matrix()");
m[i] -= ncl;
}
return m;
}

Я думаю, что это связано с показателями, потому что с меньшими числами все в порядке.
спасибо за любую помощь! я был бы так счастлив, если бы все работало так, как должно.

Джефф