using global SEXP to store XPtr gives garbage values

Question

I am building a package which solves Ordinary Differential Equations (ODEs) using the CVODE C routine (part of the SUNDIALS C library).

The package works if the user supplies a function which calculates derivatives and has the following form

#include <Rcpp.h>
using namespace Rcpp;

#include <cvode/cvode.h>               /* prototypes for CVODE fcts., consts. */
#include <nvector/nvector_serial.h>    /* serial N_Vector types, fcts., macros */

int test (realtype t, N_Vector y, N_Vector ydot, void *user_data){

  // test function
  NV_Ith_S(ydot,0) = 1*NV_Ith_S(y,0);
  NV_Ith_S(ydot,1) = 2*NV_Ith_S(y,1);
  NV_Ith_S(ydot,2) = 3*NV_Ith_S(y,2);

  return(0);

} 

typedef int (*funcPtr)(realtype t, N_Vector y, N_Vector ydot, void *user_data);

// [[Rcpp::export]]
XPtr<funcPtr> putFunPtrInXPtr() {

  // return(XPtr<funcPtr> (new funcPtr(&test)));
  XPtr<funcPtr> testptr(new funcPtr(&test), false);
  return testptr;

A function pointer, i.e. my_fun <- putFunPtrInXPtr() is formed in R and my_fun is provided to the cvode function from the package (cvode inputs my_fun as an SEXP, see code here). This works, i.e. gives the right results (see detailed instructions here). However, this requires that user has SUNDIALS installed on their system (to access cvode.h and nvector_serial.h).

I am trying to make the package so the user does not need to have SUNDIALS installed. So, the function to get derivatives (and generate function pointer) will be as below

#include <Rcpp.h>
using namespace Rcpp;

//---------------------------------------------------------------------------------
typedef NumericVector (*funcPtr1) (double t, NumericVector y, NumericVector ydot);
//---------------------------------------------------------------------------------

// [[Rcpp::export]]
NumericVector test1 (double t, NumericVector y, NumericVector ydot){

  ydot[0] = 1 * y[0];
  ydot[1] = 2 * y[1];
  ydot[2] = 3 * y[2];

  return ydot;

}

// [[Rcpp::export]]
XPtr<funcPtr1> putFunPtrInXPtr1() {

    XPtr<funcPtr1> testptr1(new funcPtr1(&test1), false);
    return testptr1;

}

On the R side, my_fun1 <- putFunPtrInXPtr1() will be run and my_fun1 will be provided to cvode_test (a test function defined in package to be able to handle derivative functions defined using NumericVector only.

In my package to convert function pointer type from XPtr<funcPtr1>to XPtr<funcPtr>, I do the following

1) A global SEXP (sexp_g) is defined outside any function 2) In cvode_test the input SEXP is assigned to sexp_g 3) Finally, a function fun_test1 is defined as follows, which converts N_Vector to NumericVector, uses the function in sexp_g to get derivatives and then puts them again in an N_Vector, i.e.

typedef int (*funcPtr_test)(double time, NumericVector y, NumericVector ydot);

SEXP sexp_g;  // declare a global SEXP

int fun_test1(realtype t, N_Vector y, N_Vector ydot, void* user_data){

  // convert y to NumericVector y1
   int y_len = NV_LENGTH_S(y);

   NumericVector y1(y_len);    // filled with zeros
   for (int i = 0; i < y_len; i++){
     y1[i] = NV_Ith_S(y,i);
   }

  // use function pointer to get the derivatives
  XPtr<funcPtr_test> xpfun(sexp_g);
  funcPtr_test fun_test = *xpfun;

  NumericVector ydot1(y1.length());
  ydot1 = fun_test(t, y1, ydot1);

  // convert ydot1 to N_Vector ydot
  // N_Vector ydot; ydot = NULL;
  ydot = N_VNew_Serial(ydot1.length());
  for (int i = 0; i<ydot1.length(); i++){
    NV_Ith_S(ydot, i) = ydot1[i];
  }


  return (0);
}

Finally, in cvode_test, this fun_test1 is used as follows

flag = CVodeInit(cvode_mem, fun_test1, T0, y0);

This also compiles but the issue is when I supply my_fun1 (i.e., a pointer to test1) to cvode_test for integration, I get garbage values back. I am not sure what is going wrong here, I have read a few articles about protecting SEXP (i.e., the one here) but I don't know how to implement it here.

Any help regarding what is wrong here, and if there is a better approach than a global SEXP variable will be helpful. Since the rhs function for CVOdeInit has to have the following signature

int RHS_function (N_Vector, N_Vector, void*);

I am finding it impossible to insert information in a function defined as above from a function defined outside the package using a different signature (using an SEXP), other than using a global variable and have been struggling with this problem for quite some time. Any help would be much appreciated!

A different approach that I have not tried is 1) declaring XPtr<funPtr_test> outside any function, 2) unwrapping SEXP to XPtr<funcPtr_test> inside the cvode_test and assigning it to the pointer declared outside the function. But I am struggling with the syntax to declare a function pointer of the type XPtr<funcPtr_test>.

The full code for cvode and cvode_test can be found here

Thanks