Fw : Re: segmentation fault using qag.h with GSL

[chadia kanaan <chadiakanaan@........>]

--- En date deÂ: Mar 12.8.08, chadia kanaan <chadiakanaan@........> a ÃcritÂ:
De: chadia kanaan <chadiakanaan@........>
Objet: Re: [Help-gsl] segmentation fault using qag.h with GSL
Ã: "Thomas Weber" <thomas.weber.mail@.........>
Date: Mardi 12 AoÃt 2008, 11h49

Dear Thomas,

I re-enclose again the error file concerning the debugging in a plain text format ( for the convenience of kind users that may help me ). You'll find also attached the code file "code.c" .
I really appreciate your help, many thanks, also for Joel :
" Putting it in MS Word adds nothing to understanding your code", surely true, 
but made me understand and admire your profound dedication to the GNU project !!

All the best of luck,
Chadia





--- En date deÂ: Mar 12.8.08, Thomas Weber <thomas.weber.mail@.........> a ÃcritÂ:
De: Thomas Weber <thomas.weber.mail@.........>
Objet: Re:
 [Help-gsl] segmentation fault using qag.h with GSL
Ã: help-gsl@gnu.org
Cc: "chadia kanaan" <chadiakanaan@........>
Date: Mardi 12 AoÃt 2008, 9h39

Am Montag,
 den 11.08.2008, 17:02 +0000 schrieb chadia kanaan:
> Dear all,
> 
> I
 am a new user of both GSL library and Ubunto 8.04 LTS.
> I have passed the compilation stage of my code successfully, however every
time I try to run it, I get a segmentation fault on the level of the integration
function qag.
> 
> I tested the  function with the example in the reference manual and it did
work.
> I upgraded from gsl 1.10-4 to GSL.11.1 and I always get the same error.
> 
> Please, if  anybody could help me, I'll be very greatfull,
> All the best,
> Chadia
> 
> PS: the details of my compilation are attached to this message.
>        In case you need the whole code, please, tell
 me

Yes, please provide the full code.

Maybe someone can figure out the problem with the short snippet you
gave, but it usually is much easier to compile the code oneself and have
a
 look with gdb. 

	Thomas






       
Envoyà avec Yahoo! Mail.
Une boite mail plus intelligente. 


      _____________________________________________________________________________ 
Envoyez avec Yahoo! Mail. Une boite mail plus intelligente http://mail.yahoo.frkanaan@port-kanaan:~/simulation nice$ gcc -Wall -W -Werror -g -o code code.c -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm  
kanaan@port-kanaan:~/simulation nice$ ./code
Segmentation fault
kanaan@port-kanaan:~/simulation nice$ gdb
GNU gdb 6.8-debian
Copyright (C) 2008 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "i486-linux-gnu".
(gdb) file co
co: No such file or directory.
(gdb) file co
code              code.c            command_file.doc  
(gdb) file code
Reading symbols from /home/kanaan/simulation nice/code...done.
(gdb) run
Starting program: /home/kanaan/simulation nice/code 

Program received signal SIGSEGV, Segmentation fault.
0xb7e54b0d in gsl_integration_qag (f=0xbfe94b60, a=0, b=100, epsabs=0, 
    epsrel=9.9999999999999995e-08, limit=50, key=6, workspace=0x80721f8, 
    result=0xbfe94b70, abserr=0xbfe94b68) at initialise.c:31
31	initialise.c: No such file or directory.
	in initialise.c
(gdb) bt
#0  0xb7e54b0d in gsl_integration_qag (f=0xbfe94b60, a=0, b=100, epsabs=0, 
    epsrel=9.9999999999999995e-08, limit=50, key=6, workspace=0x80721f8, 
    result=0xbfe94b70, abserr=0xbfe94b68) at initialise.c:31
#1  0x08049276 in p_simfunc (z=69.695184840895408, SFR_index=1) at code.c:147
#2  0x08049359 in main () at code.c:175
(gdb) b
Breakpoint 1 at 0xb7e54b0d: file initialise.c, line 31.
(gdb) up
#1  0x08049276 in p_simfunc (z=69.695184840895408, SFR_index=1) at code.c:147
147	        gsl_integration_qag (&F, zmin, zmax, 0, 1e-7, 50, 6,w, &int_p_unor, &err);
(gdb) up
#2  0x08049359 in main () at code.c:175
175	       if (z_test > p_simfunc(z_rand, SFR_index))  
(gdb) up
Initial frame selected; you cannot go up.
(gdb) 


**************************************************************************************
NOTE : the file HEADER are installed in usr/local/include
       the libraries are installed in usr/local /lib  (see below)

kanaan@port-kanaan:/usr/local/lib$ ls
libgsl.a        libgslcblas.so        libgsl.la    libgsl.so.0.12.0  
libgslcblas.la  libgslcblas.so.0.0.0  libgsl.so.0 
*****************************************************************************************
the piece of code where the compiler is stuck :
---------------------------------------------
/*define the func to be integrated*/
double inv_H(double z, void *params_ptr) 
      {
      double params = *(double*)params_ptr;
      /*printf("H_inv test passed succesfully: %g\n", params*pow(Hfunc(z), -1.0)); */
      double inv_H = params*pow(Hfunc(z), -1.0);  
      return inv_H;
      }

/*define Dl*/
double Dlfunc(double z)
       { 
        double Dl;
        double int_invH, err; /*declare the ouptput of the integration and its error*/
        double params = 1.0;
        gsl_integration_workspace* w=gsl_integration_workspace_alloc(10000);
	
	F.function = &inv_H;
	F.params = &params;
        gsl_integration_qag (&F, zmin, z, 0, 1e-7, 50, 6, w, &int_invH, &err);
        Dl=c*(1+z)*int_invH ;
        gsl_integration_workspace_free (w);/*clean workspace for the second call of Dlfunc*/
        return Dl;
       }



#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_integration.h>

#define H0 70
#define omegav 0.7
#define omegam 0.3
#define omegak 0.0
#define zmax 100
#define zmin 0
#define conv_factor 1
#define c 1     

/************************** List of FUNCTIONS *********************************************/

double Hfunc(double z)
          {  
           double H;
	   H=H0*pow((omegam*pow((1+z),3)+omegak*pow((1+z),2)+omegav),0.5);
           /*printf("hubble test passed succesfully: %g\n",H);*/
           return H;
          }
  
double SFRfunc(double z, int SFR_index)
      {  
      double SFR;
      if (SFR_index == 1)
              {
             SFR=0.32*(exp(3.4*z)/(exp(3.8*z)+45))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==2)
              {
             SFR=0.16*(exp(3.4*z)/(exp(3.4*z)+22))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==3)
              {
             SFR=0.22*(exp(3.05*z-0.4)/(exp(2.93*z)+15))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==4) 
              { 
                if ((z>=0) && (z<=1.04))
                   {
		   SFR=0.015*pow((1+z),3.28)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }
                else if ((z>=1.04) && (z<=4.48))
                   {
                   SFR=0.188*pow((1+z),-0.26)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }
                else /*if ((z>=4.48) && (z<=6))*/
                   {
                   SFR=97.72*pow((1+z),-8.0)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }                          
              }

      else if (SFR_index==5)
              {
             SFR=(0.2331*(pow((1+z),3.7))/pow(1+0.075*pow(1+z,3.7),1.84))*Hfunc(z);
              } 
      else if (SFR_index==6)
              {
            SFR=(0.0103+0.088*z)/((1+pow(z/2.4,2.8))*pow(0.3*pow((1+z),3)+0.7,0.5))*Hfunc(z);
              } 
      else  if (SFR_index==7)
              {
             SFR=(0.014+0.11*z)/((1+pow(z/1.4,2.2))*pow(0.3*pow((1+z),3)+0.7,0.5))*Hfunc(z);
              } 
      else 
              { 
                if ((z>=0) && (z<=0.993))
                   {
                   SFR=0.02*pow(1+z,3.30)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                else if ((z>=0.993) && (z<=3.80))
                   {
                   SFR=0.19*pow(1+z,-0.0549)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                else /*if ((z>=3.80) && (z<=10))*/
                   {
                   SFR=223.9*pow(1+z,-4.46)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                /*else break*/                        
              }
              printf("SFR test passed succesfully: %g\n",SFR);       
   return SFR;
      } 


double inv_H(double z, void *params_ptr) 
      {
      double params = *(double*)params_ptr;
      /*printf("H_inv test passed succesfully: %g\n", params*pow(Hfunc(z), -1.0)); */
      double inv_H = params*pow(Hfunc(z), -1.0);  
      return inv_H;
      }


double Dlfunc(double z)
       { 
        double Dl;
        double int_invH, err; /*declare the ouptput of the integration and its error*/
        double params = 1.0;
        gsl_integration_workspace* w=gsl_integration_workspace_alloc(10000);
	/*gsl_function F = {&H_inv,0};*/ 
        gsl_function F;
	F.function = &inv_H;
	F.params = &params;
        /*printf("integration_before_test passed succesfully")*/
        gsl_integration_qag (&F, zmin, z, 0, 1e-7, 50, 6, w, &int_invH, &err);
        Dl=c*(1+z)*int_invH ;
        printf("Dl test passed : %g\n", Dl);
        gsl_integration_workspace_free (w);/*clean workspace for the second call of Dlfunc*/
        return Dl;
       }


double dvdzfunc(double z)
       { 
        double dvdz;
        dvdz=4*M_PI*c*pow(Dlfunc(z),2)/(pow(1+z,2)*Hfunc(z));
        return dvdz;
       }


double p_unor( double z, void *SFR_index_ptr) 
{
  int SFR_index = *(int*)SFR_index_ptr;
  return conv_factor*SFRfunc(z, SFR_index)*dvdzfunc(z)*pow(1+z,-1);
}


double p_simfunc(double z, int SFR_index)
       {
        double p_sim ;
        double int_p_unor, err;
        /*normalize p_unor*/
        gsl_integration_workspace*w=gsl_integration_workspace_alloc(10000);
        gsl_integration_workspace_free(w);
	gsl_function F = { &p_unor, &SFR_index };
        gsl_integration_qag (&F, zmin, zmax, 0, 1e-7, 50, 6,w, &int_p_unor, &err);
        p_sim=p_unor(z,&SFR_index)/int_p_unor;
        gsl_integration_workspace_free (w);
        return p_sim;
       }
       
/************************************* MAIN ***************************************************/
int main (void)
{
 /********************************* simulated REDSHIFT *************************************/
   int SFR_index=1;       
   int i,i_sim=1; 
   int i_max=10000;
   double z_rand, z_test, z_accept;
   printf("program started");
   const gsl_rng_type * T;  
   gsl_rng * r;             
   gsl_rng_env_setup();     
   T = gsl_rng_ranlxd2;     
   r = gsl_rng_alloc (T);  
   gsl_vector *u=gsl_vector_alloc (10000);  
   gsl_vector *g=gsl_vector_alloc (10000);                                                     
   printf("generator declared");
   
   for(i=0;i<i_max;i++)
     {     /*use Rejection method*/
       z_test=gsl_rng_uniform(r);  
       z_rand=gsl_rng_uniform(r)*zmax;  
       if (z_test > p_simfunc(z_rand, SFR_index))  
       	 break;
       printf("rejection 1 test passed succesfully");
       z_accept = z_rand;  
       gsl_vector_set(g,i,z_rand); 
       gsl_vector_set(u,i_sim,z_accept); 
       i_sim = i_sim +1;
     }
   for(i=0;i<i_max;i++) /*display simulated redshift distribution before  rejection*/
     {
       printf("%f \n", gsl_vector_get(g,i));
     }
   for(i=0;i<i_sim;i++) /*display simulated redshift distribution after rejection*/
     {
       printf("%f \n", gsl_vector_get(u,i));
     }
     
   gsl_rng_free(r);  /*cleanup the generator workspace, and the input vectors*/
   gsl_vector_free(u);
   gsl_vector_free(g);
   return 0 ;
}