Velocity Reviews > C++ > Segmentation fault using NUMREC, MUST BE EASY!!!

# Segmentation fault using NUMREC, MUST BE EASY!!!

wollvieh
Guest
Posts: n/a

 08-27-2005
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[i][j];
cout <<"a["<<i<<"]["<<j<<"] = "<<a[i][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[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_vector(vv,1,n);
}
#undef TINY
#undef NRANSI

Victor Bazarov
Guest
Posts: n/a

 08-28-2005
wollvieh wrote:
> 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[i][j];
> [...]

V

Frank Chang
Guest
Posts: n/a

 08-29-2005
wollvieh, You can save your time by using the boost template class
Also, the boost numerical linear algebra packages are more generic than
numerical recipes.

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are Off Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post lurch132002@yahoo.com C Programming 4 05-01-2007 01:52 AM jgarber Python 0 08-16-2006 06:20 PM jtagpgmr C Programming 11 06-27-2006 06:59 PM Keith Lee Firefox 3 04-29-2006 05:45 PM cesco C++ 1 03-03-2006 03:44 PM