simple c++ help required!

AGD

AGD

Soldato
Joined
23 Nov 2007
Posts
5,048
Hi,

I've found a code to solve linear equations (Ax = b - all matrices/vectors) and need to get it running asap. I can understand all the bits inside the function but am having some trouble with the arguments. Am I correct in thinking the function takes the addresses of x and b as arguments? I don't know what the ** before the a means though.

Also can I just change the

#define eps 1e-10

to

const double eps?

Code:
#include <math.h>

#define EPS 1e-10

int gelimd2(double **a,double *b,double *x, int n)
{
    double tmp,pvt,*t,**aa,*bb;
    int i,j,k,itmp,retval;

    retval = 0;
    aa = new double *[n];
    bb = new double [n];
    for (i=0;i<n;i++) {
        aa[i] = new double [n];
        bb[i] = b[i];
        for (j=0;j<n;j++) {
            aa[i][j] = a[i][j];
        }
    }
    for (i=0;i<n;i++) {             // outer loop on rows
        pvt = aa[i][i];              // get pivot value
        if (fabs(pvt) < EPS) {
            for (j=i+1;j<n;j++) {
                if(fabs(pvt = aa[j][i]) >= EPS) break;
            }
            if (fabs(pvt) < EPS) {
                retval = 1;
                goto _100;     // pull the plug!
            }
            t=aa[j];                 // swap matrix rows...
            aa[j]=aa[i];
            aa[i]=t;
            tmp=bb[j];               // ...and result vector
            bb[j]=bb[i];
            bb[i]=tmp;        
        }
// (virtual) Gaussian elimination of column
        for (k=i+1;k<n;k++) {       // alt: for (k=n-1;k>i;k--)
            tmp = aa[k][i]/pvt;
            for (j=i+1;j<n;j++) {   // alt: for (j=n-1;j>i;j--)
                aa[k][j] -= tmp*aa[i][j];
            }
            bb[k] -= tmp*bb[i];
        }
	}
// Do back substitution
	for (i=n-1;i>=0;i--) {
        x[i]=bb[i];
		for (j=n-1;j>i;j--) {
            x[i] -= aa[i][j]*x[j];
		}
        x[i] /= aa[i][i];
	}
// Deallocate memory
_100:
    for (i=0;i<n;i++) {
        delete [] aa[i];
    }
    delete [] aa;
    delete [] bb;
    return retval;
}

Thanks.
 
Back
Top Bottom