Using params from gsl function in integration routine

[CAF123]

I'm trying to pass through some parameters of a function to the gsl integration routine but my code is currently not returning correct values. I attach a version of my code using dummy example functions and names.

struct myStruct_t {
           double a;
 };

double func(double z, void* params)  {
        myStruct_t &myStruct = *(myStruct_t *)(params);

        double res = pow(z*myStruct.a,2);

       return res;
}
     
int main() {
    double result, error;
    myStruct_t myStruct;
    gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);

    gsl_function F;
    F.params = (void*)(&myStruct);
    F.function = func;

    for(int i = 0; i < 10; i++) {
   
        double b = pow(2,i);
        double a = 2*b;

        myStruct.a = a;
        gsl_integration_qags(&F, 0, 1, 1e-3, 1e-3, 1000, w, &result, &error);
        std::cout << b << result << std::endl;
    }

}

I believe the issue is with the myStruct.a parameter in ``func and somehow it is not being passed through the integration routine properly. The code does not return an error but yields values ~10^(129) suggesting something is not quite right. I have tested this by simply hardcoding ten values of a one at a time in func instead of using a for loop and seeing if the code returns sensible results which it does so. Do you see immediately any error in what I have written? Thanks in advance!

 

[jedishrfu]

Use a debugger and step through your code. The effort to learn how to use the debugger will definitely payoff in the long run.

 

[Svein]

myStruct_t &myStruct = *(myStruct_t *)(params);

I may be wrong (I haven't analysed C code since 2010), but I think the quoted part is subtly wrong. Going from right to left:

• (myStruct_t *)(params) casts the typeless params to be a pointer to a myStruct_t
• *(myStruct_t *)(params) dereferences that pointer (in other words it fetches the value that params is pointing to)

On the left side:

• myStruct_t &myStruct is a pointer

So what you are doing is assigning a value (the contents of a struct) to a pointer

My standard advice to all C programmers: Run the code through lint repeatedly until no more comments appear.

 

[Filip Larsen]

myStruct_t &myStruct = *(myStruct_t *)(params);
[...]
On the left side:

• myStruct_t &myStruct is a pointer

The syntax on the left side of the line you quoted means that myStruct is declared as a myStruct_t reference, not a pointer, and in the given context it does appears to be correct.

I would also recommend using a debugger to verify if the params parameter really is pointing to the myStruct instance or not. If it doesn't then I guess it requires some gsl specific knowledge or documentation to figure out why not.

 

[CAF123]

Dear all, Thanks for your comments! Actually I am using the lldb debugger inbuilt into xcode on a macos. The program builds but there is 'no debug session' initiated in the IDE interface.

I had looked here at the first example
https://www.gnu.org/software/gsl/doc/html/integration.html

It's not immediately helpful because alpha is declared as a double dereferenced from pointer but then in the actual int main alpha is just set as 1.0.

Maybe @DrClaude can help, I recall I had discussed gsl with him/her before.

 

[jedishrfu]

With an IDE debugger, you run the code in debug mode. However, prior to running in debug mode, you must set a breakpoint on the line of interest or the line before.

 

[CAF123]

With an IDE debugger, you run the code in debug mode. However, prior to running in debug mode, you must set a breakpoint on the line of interest or the line before.

Thanks, I see, putting the breakpoint after the return res in double func in my actual code says that the res is returning a value of around -1e+126. But I haven't declared any value for

 a

yet (that comes later in the for loop) so how can it resume this value? Is it just like some temporary memory value or something?

 

[jedishrfu]

We think your problem is in how you're referencing the mystruct variable so that means stepping into the function and looking at what it thinks the arguments are. Sometimes these types of errors are more complex than you can imagine.

I had an issue across the calling boundary where data was getting corrupted. It turned out the library code I was calling was compiled with different flags from my code where the library code expected argument data to be placed at 8 byte boundaries instead of 4 byte.

 

[Svein]

The syntax on the left side of the line you quoted means that myStruct is declared as a myStruct_t reference, not a pointer, and in the given context it does appears to be correct.

OK. I am no expert on C++ (on the other hand I have been programming in C for 25years), but I still do not see why you would want to do it this way. In C, the obvious way would be:

double func(double z, void* params)  {
        myStruct_t myStruct = *(myStruct_t *)(params);

        double res = pow(z*myStruct.a,2);

       return res;
}

The difference being that you now create a local myStruct and copy the given values to it.

BTW - this video () lists void * as one of the 10 "forbidden" C++ constructs.

 

[pbuk]

If you look at this example: https://www.gnu.org/software/gsl/doc/html/roots.html#c.gsl_function I think you should be using

    gsl_function F;
    F.params = &myStruct;
    F.function = &func;

This example also shows how to destructure the params in func using ->.

 

[DrClaude]

params is being passed as a pointer of type void and you simply want to recast it into a pointer of type myStruct_t. The elements of the struct can then be accessed using -> instead of .

struct myStruct_t {
           double a;
};

double func(double z, void* params)  {
       // myStruct_t &myStruct = *(myStruct_t *)(params);
        myStruct_t *myStruct = (myStruct_t *)(params);
        double res = pow(z*myStruct->a,2);

       return res;
}

Maybe @DrClaude can help, I recall I had discussed gsl with him/her before.

I use GSL in my teaching, so don't hesitate to tag me when you have such questions.

 

[Filip Larsen]

I am no expert on C++ (on the other hand I have been programming in C for 25years), but I still do not see why you would want to do it this way.

It is mostly a matter of style and convention. In order not to derail this thread let me just refer to this Pointer vs Reference in C++ page for details. In context of this thread my point was that both approaches will work equally well inside the func function, even so much as to produce the same compiled code.

Granted, if one is using a C++ library that has a C-flavored API and examples, then it certainly may be a good idea to follow the "library style" when in doubt.