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]);
}
#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]);
}