help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Segfault in gsl_integration_qag


From: Thomas Markovich
Subject: [Help-gsl] Segfault in gsl_integration_qag
Date: Sun, 24 May 2009 09:00:59 -0500

Hi,

I have been racking my brain over this segmentation fault for a few days now and tracked it down to gsl_integration_qag. The strange thing is that the segmentation fault only happens after n and m = 37. I did try increasing the size of my work space but that didn't seem to do much of anything at all. Any suggestions? I'm sure the error is something thats obvious to someone more experienced with the gsl.

I have included the code below:

        int N = 50; //Starting size
        double sum = 0.0;
        int j = 0;
        int sizeN = 2*N+1; //Size
        mat HN(sizeN,sizeN); //Defines the hamiltonian matrix
gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);
        double result1,result2, error;
        gsl_function F1;
        struct my_f_params alpha = {2,2};               
        F1.function = &f1;
        F1.params = α
        for(int n = 0; n <= sizeN; n++)
        for(int m = 0; m <= sizeN; m++)
        {
                        alpha.a = n;
                        alpha.b = m;
gsl_integration_qag (&F1, 0, 6, 0, 1e-4, 1000, GSL_INTEG_GAUSS61,w, &result1, &error); gsl_integration_qag (&F1, -6, 0, 0, 1e-4, 1000, GSL_INTEG_GAUSS61,w, &result2, &error);
                        cout << m << "  " << n << endl;
                        HN(n,m) = result1 + result2;
        }

with f1 being

double f1 (double x, void * p) { //Deriviative part of the integral
        struct my_f_params * params = (struct my_f_params *)p;
        int n = (params->a);
        int m = (params->b);
double f = exp(-(x*x))*hermitecalculate(m,x)*( (-1.0)*((4*n*n - 4*n)*hermitecalculate(n-2,x) - 4*n*hermitecalculate(n-1,x)*x - hermitecalculate(n,x) + hermitecalculate(n,x)*x*x)); f += exp(-x*x)*hermitecalculate(n,x)*hermitecalculate(m,x)*(pow(x, 6.0) + 4*pow(x,4) + x*x - 2); return f/ (sqrt (gsl_sf_fact (n)*sqrt(pi)*pow(2.0,n))*sqrt(gsl_sf_fact(m)*sqrt(pi)*pow(2.0,m)));
}

and my_f_params defined as

struct my_f_params {int a; int b;};




reply via email to

[Prev in Thread] Current Thread [Next in Thread]