Octave/C++ matrix Inv() comparison

[ionone]

Hi

i'm trying to convert an Octave .m file to C++. So far i quite succeded
except for a matrix inverse function : i don't have the same numbers in
Octave and in my C++ program. i checked, double checked all the values and i
cannot find a reason why the two values differ.

I can't post here a SSCCE because i couldn't find one but i can post the
function i use (at the end of this message)

the kind of values i get is for example in Octave :3.8070e+010 and in C++:
38121149284.106712

or -1.4971e+011 in Octave and -149962686456.46307 in C++

so there is a difference and i messes all up in my function...

please please if you don't know give me directives so i can progress because
two days i'm on this single issue and i'm really blocked.

thanks

Jeff



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, &amp;d);
  for (j=0; j&lt;n; j++) 
  {
    d *= a[j][j];
    for (i=0; i&lt;n; i++) 
        col[i]=0.0;
    col[j]=1.0;
    lubksb(a,n,indx,col);
    for (i=0; i&lt;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&lt;n;i++) {
    big=0.0;
    for (j=0; j&lt;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&lt;j; i++) {
      sum=aa[i][j];
      for (k=0; k&lt;i; k++) 
          sum -= aa[i][k]*aa[k][j];
      aa[i][j]=sum;
    }
    big=0.0;
    for (i=j; i&lt;n; i++) 
    {
      sum=aa[i][j];
      for (k=0; k&lt;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 &lt; 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&lt;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;
}



--
View this message in context: 
http://octave.1599824.n4.nabble.com/Octave-C-matrix-Inv-comparison-tp4655291.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.