Segmentation fault using NUMREC, MUST BE EASY!!!

Discussion in 'C++' started by wollvieh, Aug 27, 2005.

  1. wollvieh

    wollvieh Guest

    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
     
    wollvieh, Aug 27, 2005
    #1
    1. Advertising

  2. 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[j];
    > [...]


    V
     
    Victor Bazarov, Aug 28, 2005
    #2
    1. Advertising

  3. wollvieh

    Frank Chang Guest

    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.
     
    Frank Chang, Aug 29, 2005
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. jgarber
    Replies:
    0
    Views:
    722
    jgarber
    Aug 16, 2006
  2. jtagpgmr

    Segmentation Fault using the srand function

    jtagpgmr, Jun 27, 2006, in forum: C Programming
    Replies:
    11
    Views:
    981
    Default User
    Jun 27, 2006
  3. cesco
    Replies:
    1
    Views:
    710
  4. Replies:
    4
    Views:
    343
    Ian Collins
    May 1, 2007
  5. Replies:
    5
    Views:
    542
    santosh
    Aug 26, 2007
Loading...

Share This Page