[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Octave/C++ matrix Inv() comparison
From: |
ionone |
Subject: |
Octave/C++ matrix Inv() comparison |
Date: |
Sun, 7 Jul 2013 08:36:18 -0700 (PDT) |
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, &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;
}
--
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.