Need help passing array to ludcmp


apwinter
10-27-2008, 02:58 PM
Hi, I am new to C code and I am trying to use ludcmp and lubskb. I am able to get it to compile, but the vector it returns (b) is not correct. Can someone look at the code and let me know what I am doing incorrect? I have spent a lot of time on this and can not figure it out. The routines ludcmp and lubskb are straight from NR with the exception I am using double instead of float. Could there also be a correction to the NR routines that I have missed from the forums?

#include <math.h>
#include "nrutil.h"
#include "stdio.h"
#define TINY 1.0e-20

void ludcmp(double **a, int n, int *indx, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;
vv=dvector(1,n);

*d=1.0;
for (i=1; i<=n; i++) {
big=0.0;
for (j=1; j<=n; j++)
if ((temp=fabs(a[i][j]))>big) big=temp;
if (big==0.0) nrerror("Singular matrix in routine ludcmp");
vv[i]=1.0/big;
}
for (j=1;j<=n;j++) {
for (i=1;i<j;i++) {
sum=a[i][j];
for(k=1;k<i;k++) sum -= a[i][k] * a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<=n;i++) {
sum=a[i][j];
for (k=1;k<j;k++)
sum -= a[i][k] *a [k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >=big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=1;k<=n;k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}

indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != n) {
dum=1.0/(a[j][j]);
for (i=j+1;i<=n;i++) a[i][j] *= dum;
}
}
free_dvector(vv,1,n);
}

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

for (i=1;i<=n;i++) {
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=1;
b[i]=sum;
}
for (i=n;i>=1;i--) {
sum=b[i];
for (j=i+1;j<=n;j++) sum -= a[i][j] * b[j];
b[i]=sum/a[i][j];
}
}

main()
{
double d, **a, *b;
int n=3,*indx;

indx=ivector(1,n);
a=dmatrix(1,n,1,n);
b=dvector(1,n);

a[1][1] = 2.;
a[1][2] = 1.;
a[1][3] = -2.;
a[2][1] = 1.;
a[2][2] = -5.;
a[2][3] = 3.;
a[3][1] = 3.;
a[3][2] = -3.;
a[3][3] = -2.;
b[1]=8.;
b[2]=6.;
b[3]=10.;

ludcmp(a,n,indx,&d);

lubksb(a,n,indx,b);

printf("%f ",b[1]);
printf("%f ",b[2]);
printf("%f \n",b[3]);

}

davekw7x
10-27-2008, 03:29 PM
...the vector it returns (b) is not correct.
You have a couple of transcription errors in your lubksb() function. See Footnote.

Your line 77

ii = i; /* It was ii = 1; */


Your line 84

b[i] = sum / a[i][i];/* It was b[i] = sum / a[i][j]; */


Answer:

b[1] = 5.200000, b[2] = 0.800000, b[3] = 1.600000


Regards,

Dave


Footnote from the text:

...experience has shown that virtually all reported bugs in such cases are typing errors!

apwinter
10-27-2008, 03:49 PM
Thank you so much. I can not believe I overlooked the typos. My mind is at ease now. I can move on to something else in my research. Thank you for responding so quickly.