2015-09-04 175 views
3

我試圖從數值配方(代碼取自here,第164頁)實現多維集成例程。使用包裝來傳遞被積函數到void* params是從here(進一步參考那裏)。我class聲明在NRquad3d_Const_qags_2.C如下調用gsl_integration時出現分段錯誤

#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   = % .18f\n", result); 
    printf ("exact result = % .18f\n", expected); 
    printf ("actual error = % .18f\n", result - expected); 

    return 0; 
} 

代碼編譯好了,但運行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已經幫我消除在w1w2地址的衝突,但段錯誤仍然存​​在。由於錯誤似乎發生在第二次調用qagsNRf2的調用中,我認爲必須有一個釋放的指針,我試圖引用,但我似乎無法找到它。誰能幫忙?

回答

0

運行valgrind給出了形式的錯誤:

==26286== Invalid read of size 8 
==26286== at 0x4ED0239: gsl_integration_workspace_free (in /usr/lib64/libgsl.so.0.16.0) 
==26286== by 0x400A54: main (crap_test.cpp:32) 
==26286== Address 0x6053090 is 80 bytes inside a block of size 88 free'd 
==26286== at 0x4C29D4E: free (in /usr/lib64/valgrind/vgpreload_memcheck-amd64-linux.so) 
==26286== by 0x400C39: NRf2::ret_int_2(double) (NRquad3d_Const_qags_2.C:46) 
==26286== by 0x400C76: NRf2::ret_int_2_wrapper(double, void*) (NRquad3d_Const_qags_2.C:53) 
==26286== by 0x4ECBBA3: gsl_integration_qk (in /usr/lib64/libgsl.so.0.16.0) 
==26286== by 0x4ECB97A: gsl_integration_qk21 (in /usr/lib64/libgsl.so.0.16.0) 
==26286== by 0x4ECD4C6: ??? (in /usr/lib64/libgsl.so.0.16.0) 
==26286== by 0x400D51: output::integrate(double, double, double, double, double) (NRquad3d_Const_qags_2.C:73) 
==26286== by 0x4009CB: main (crap_test.cpp:25) 

這暗示我trying to access freed memory。在我的腳本中的唯一的地方我嘗試這樣做是在線路

gsl_integration_workspace_free (w1); 
. 
. 
. 
. 
gsl_integration_workspace_free (w2); 
.C文件

。然後我(重新)將這些行移動到我的.cpp文件的末尾,它可以工作。

愚蠢的錯誤,但我會回答這個,無論如何有人有這個確切的(或類似的!)錯誤,需要一個解決方案。

+1

如果你真的進入C++,使用unique_ptr與自定義析構函數不再有這個問題了! –