Search code examples
c++classshared-librariesgsl

Calling GSL function inside a class in a shared library


I'm trying make a shared library in c++ implementing tools for Fermi gases. I'm using the GSL library to solve a function numerically and my code runs without a problem without when running as a script, but when trying to convert it to a shared library and classes I encounter problems.

I've seen similar questions: Q1 Q2 Q3

I'm fairly new to c++-programming and cannot seem to adapt the different answers to my problem. Probably since I do not quite understand the answers.

My code is:

/* Define structure for the GSL-function: chempot_integrand */
struct chempot_integrand_params { double mu; double T; };

double
ChemicalPotential::chempot_integrand (double x, void * params){
    /* Computes the integrand for the integral used to obtain the chemical potential.
     *
     * This is a GSL-function, which are integrated using gsl_integration_qag.
     */

    // Get input parameters.
    struct chempot_integrand_params * p = (struct chempot_integrand_params *) params;
    double mu = p->mu;
    double T = p->T;

    // Initiate output parameters for GSL-function.
    gsl_sf_result_e10 result;
    int status = gsl_sf_exp_e10_e( ( gsl_pow_2(x) - mu ) / T , &result );

    if (status != GSL_SUCCESS){
        printf ("Fault in calculating exponential function.");
    }

    // Return (double) integrand.
    return (gsl_pow_2(x) / ( 1 + result.val * gsl_sf_pow_int(10,result.e10) ));
}

/* Define structure for the GSL-function: chempot_integration */
struct chempot_integral_params { double T; };

double
ChemicalPotential::chempot_integration (double mu, double T){
    /* Computes the integral used to obtain the chemical potential using the integrand: chempot_integrand.
    */

    // Set input parameters for the integrand: chempot_integrand.
    struct chempot_integrand_params params_integrand = { mu, T };

    // Initiate the numerical integration.
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); // Allocate memory for the numerical integration. Can be made larger if neccessary, REMEMBER to change it in the function call: gsl_integration_qag as well.
    double result, error;
    gsl_function F;
    F.function = &ChemicalPotential::chempot_integrand;
    F.params = &params_integrand;

    // Upper limit for integration
    double TOL = 1e-9;
    double upp_lim = - T * gsl_sf_log(TOL) + 10;

    gsl_integration_qag (&F, 0, upp_lim, 1e-12, 1e-12, 1000, 6, w, &result, &error);

    // Free memory used for the integration.
    gsl_integration_workspace_free (w);

    return result;
}

and when compiling I get the error

error: cannot convert ‘double (Fermi_Gas::ChemicalPotential::*)(double, void*)’ to ‘double (*)(double, void*)’ 

in line

F.function = &ChemicalPotential::chempot_integrand;

Solution

  • It is indeed interesting that people ask this over and over again. One reason may be that the proposed solutions are not easy to understand. I for one had problems understanding and implementing them. (the solutions did not work out of the box for me, as you might expect.)

    With the help of tlamadon I just figured out a solution that may be helpful here as well. Let's see what you guys think.

    So just to recap, the problem is that you have a class that contains a member function on which you want to operate with something from the GSL library. Our example is useful if the GSL interface requires a

    gsl_function F;
    

    see here for a definition.

    So here is the example class:

    class MyClass {
    
        private:
            gsl_f_pars *p;  // not necessary to have as member
    
        public: 
            double obj(double x, void * pars);  // objective fun
            double GetSolution( void );  
            void setPars( gsl_f_pars * xp ) { p = xp; };
            double getC( void )  ;  // helper fun
    
    };
    

    The objective of this exercise is to be able to

    1. initiate MyClass test,
    2. supply it with a paramter struct (or write a corresponding constructor), and
    3. call test.GetSolution() on it, which should return whatever the GSL function was used for (the minimum of obj, a root, the integral or whatever)

    The trick is now to put have an element in the parameter struct gsl_f_pars which is a pointer to MyClass. Here's the struct:

    struct gsl_f_pars {
        double a;
        double b;
        double c;
        MyClass * pt_MyClass;
    };
    

    The final piece is to provide a wrapper that will be called inside MyClass::GetSolution() (the wrapper is a stand in for the member function MyClass::obj, which we cannot just point to with &obj inside the class). This wrapper will take the parameter struct, dereference pt_MyClass and evaluate pt_MyClass's member obj:

    // Wrapper that points to member function
    // Trick: MyClass is an element of the gsl_f_pars struct
    // so we can tease the value of the objective function out
    // of there.
    double gslClassWrapper(double x, void * pp) {
        gsl_f_pars *p = (gsl_f_pars *)pp;
        return p->pt_MyClass->obj(x,p);
    }
    

    The full example is a bit too long to post here, so I put up a gist. It's a header file and a cpp file, it should be working wherever you have GSL. Compile and run with

    g++ MyClass.cpp -lgsl -o test
    ./test