Segmentation fault using NUMREC, MUST BE EASY!!!

[wollvieh]

Hello,

this one must be easy. I need to solve a set of linear equation and
would like to use the LU decomposition ludcmp from numerical recipes. I
run into a segmentation fault and don't know why. I'm sure this is very
easy to understand for a programmer with some experience.

Here is the code:

#include <iostream>
#include <fstream>
#include <math.h>
#include "nr.h" // declare ludcmp
#include <cstdlib>
#include "nrutil.h" // declare matrix
using namespace std;

int main ()
{
float **a,*b, d;
int n, *indx;
n=2;
a=matrix(1,n,1,n);
for (int i=1;i<=n;++i)
{
for (int j=1;j<=n;++j)
{
cin >> a[j];
cout <<"a["<<i<<"]["<<j<<"] = "<<a[j]<<endl;
}
}
cout<<"Hello?";
ludcmp(a,n+1,indx,&d);
}

Reading the matrix works. But as soon as I add the last line,
ludcmp(a,n+1,indx,&d);, into the code, it returns a segmentations
fault.

Can anyone help?

The ludcmp code is, see NUMREC,

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

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

vv=vector(1,n);
*d=1.0;
for (i=1;i<=n;i++) {
big=0.0;
for (j=1;j<=n;j++)
if ((temp=fabs(a[j])) > big) big=temp;
if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
vv=1.0/big;
}
for (j=1;j<=n;j++) {
for (i=1;i<j;i++) {
sum=a[j];
for (k=1;k<i;k++) sum -= a[k]*a[k][j];
a[j]=sum;
}
big=0.0;
for (i=j;i<=n;i++) {
sum=a[j];
for (k=1;k<j;k++)
sum -= a[k]*a[k][j];
a[j]=sum;
if ( (dum=vv*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[j] *= dum;
}
}
free_vector(vv,1,n);
}
#undef TINY
#undef NRANSI

 

[Victor Bazarov]

this one must be easy. I need to solve a set of linear equation and
would like to use the LU decomposition ludcmp from numerical recipes.

And what language does 'numerical recipes' use? Not Pascal?

I run into a segmentation fault and don't know why. I'm sure this is
very easy to understand for a programmer with some experience.

Just one comment that may help you understand the problem and solve it.
Arrays in C++ are indexed from 0 (not from 1 like in Fortran). That
means that an array of 'n' elements has element indices spanning from
0 to (n-1). I see you have plenty of loops with counters going from 1
to n, which then is used to index within an array. Since I do not see
how those arrays are allocated, it's impossible to conclude with
certainty that indexing _beyond_ the [usual] end is the cause, but it
should be looked at.

Here is the code:

#include <iostream>
#include <fstream>
#include <math.h>
#include "nr.h" // declare ludcmp
#include <cstdlib>
#include "nrutil.h" // declare matrix
using namespace std;

int main ()
{
float **a,*b, d;
int n, *indx;
n=2;
a=matrix(1,n,1,n);
for (int i=1;i<=n;++i)
{
for (int j=1;j<=n;++j)

In C++ we _usually_ write

for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)

{
cin >> a[j];
[...]

V

 

[Frank Chang]

wollvieh, You can save your time by using the boost template class
<boost/numeric/ublas/lu.hpp> . Please see the link at www.boost.org.
Also, the boost numerical linear algebra packages are more generic than
numerical recipes.