调用gsl_integration时出现分段错误

我正试图从数值食谱(代码取自这里，第164页)实现多维积分例程。从这里开始使用包装器将被积函数传递给void* params(在那里有进一步的参考)。我的class声明在NRquad3d_Const_qags_2.C中如下所示

#include <gsl/gsl_integration.h>
gsl_integration_workspace * w1 = gsl_integration_workspace_alloc (100000);
gsl_integration_workspace * w2 = gsl_integration_workspace_alloc (100000);
double func3d(double x, double y, double z);
class NRf3 {
public:
  double xsav, ysav;
  double ret_int(double z)
  {  
    return func3d(xsav, ysav, z);
  }
  static double ret_int_wrapper(double z, void *params)
  {
   return static_cast<NRf3*>(params)->ret_int(z); 
  }
};
//integrates over the z-variable
class NRf2 {
private:
  double z1, z2;
  gsl_function F2; 
public:
  NRf3 f3;
  double resultf2, errorf2;
  NRf2 (double zz1, double zz2): z1(zz1), z2(zz2) {}
  double ret_int_2(double y) 
  {
    f3.ysav = y;        
    F2.function = &(f3.ret_int_wrapper);
    F2.params = &f3;
    gsl_integration_qags(&F2, z1, z2, 0, 1e-7, 10, w1, &resultf2, &errorf2);
    gsl_integration_workspace_free (w1);
    return resultf2; 
    }
  static double ret_int_2_wrapper(double y, void *params)
  {
   return static_cast<NRf2*>(params)->ret_int_2(y); 
  }
};
//integrates over the y-variable ==> gives function of x
class output {
private:  
  gsl_function F;
public:
  double result, error;
  double integrate(double x, double y1, double y2, double z1, double z2){
    NRf2 f2(z1, z2);
    f2.f3.xsav = x;
    F.function = &(f2.ret_int_2_wrapper);
    F.params = &f2;
    gsl_integration_qags(&F, y1, y2, 0, 1e-7, 10, w2, &result, &error);
    gsl_integration_workspace_free (w2);
    return result;
  }
};

并按如下方式调用积分器:

#include <stdio.h>
#include </users/.../NRquad3d_Const_qags_2.C>
double func3d(double x, double y, double z)
  {
    double f = x*y*z;
    return f;
  }

int main (void)
{
  double result; 
  double error;
  double expected = -4.0;
  output out ; 
  result = out.integrate(1.0, -1.0, 1.0, -1.0, 1.0);
  printf ("result          = % .18fn", result);
  printf ("exact result    = % .18fn", expected); 
  printf ("actual error    = % .18fn", result - expected);
  return 0;
}

代码编译ok，但运行a.out给我一个段错误。从gdb回溯得到

#0  0x00007ffff7a4b430 in ?? () from /usr/lib64/libgsl.so.0
#1  0x0000000000400bed in NRf2::ret_int_2 (this=0x7fffffffdbd0, y=-0.97390652851717174) at /users/.../NRquad3d_Const_qags_2.C:45
#2  0x0000000000400c39 in NRf2::ret_int_2_wrapper (y=-0.97390652851717174, params=0x7fffffffdbd0) at /users/.../NRquad3d_Const_qags_2.C:53
#3  0x00007ffff7a49c6e in gsl_integration_qk () from /usr/lib64/libgsl.so.0
#4  0x00007ffff7a4997b in gsl_integration_qk21 () from /usr/lib64/libgsl.so.0
#5  0x00007ffff7a4b4c7 in ?? () from /usr/lib64/libgsl.so.0
#6  0x0000000000400d14 in output::integrate (this=0x7fffffffdc30, x=1, y1=-1, y2=1, z1=-1, z2=1) at /users/.../NRquad3d_Const_qags_2.C:73
#7  0x00000000004009ac in main () at crap_test.cpp:22
(gdb) frame 6
#6  0x0000000000400d14 in output::integrate (this=0x7fffffffdc30, x=1, y1=-1, y2=1, z1=-1, z2=1) at /users/.../NRquad3d_Const_qags_2.C:73
73          gsl_integration_qags(&F, y1, y2, 0, 1e-7, 10, w2, &result, &error);
(gdb) frame 1
#1  0x0000000000400bed in NRf2::ret_int_2 (this=0x7fffffffdbd0, y=-0.97390652851717174) at /users/.../NRquad3d_Const_qags_2.C:45
45          gsl_integration_qags(&F2, z1, z2, 0, 1e-7, 10, w1, &resultf2, &errorf2);

运行gdb已经帮助我消除了w1和w2地址的冲突，但是段故障仍然存在。由于错误似乎发生在NRf2中第二次调用qags的调用中，我认为必须有一个已释放的指针，我试图引用，但我似乎找不到它。有人能帮忙吗?