Using Besselk function in C MexFunction

Question

I would like to implement the Matern correlation function in C mexFunction, which requires the computation for the modified Besselk function of the second kind.

In MATLAB, one can use the function besselk. However, there is no such equivalent in any C library (am I right?). I knew that the boost library (a C++ library) provides the access to the modified Besselk function of the second kind, see https://en.cppreference.com/w/cpp/experimental/special_math, but I have trouble making it work in my C mexFunction with MATLAB 2018a on my mac as well as a linux system. BTW, I do not want to call the MATLAB function besselk inside of C mex code.

Any suggestions will be appreciated. Below is a minimal example to construct the matern correlation function with the C mex code.

#include "mex.h"
#include "matrix.h"
#include <math.h>
//#include <boost/math/special_functions.hpp>
double matern(double h, double nu)
{
    double tau = sqrt(2.0*nu) * h;
    double c = 0.0;
    if(h>0.0){
      c = (pow(tau, nu) * pow(2.0, 1.0-nu) / tgamma(nu)) * boost::math::cyl_bessel_k(nu, tau);
    }else{
      c = 1.0;
    }
    return c;
  }
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    if(nrhs!=2){
        mexErrMsgTxt("Two Inputs are required.");
    }
    double h = mxGetScalar(prhs[0]);
    double nu = mxGetScalar(prhs[1]);
    double corr = matern(h, nu);
    mexPrint("matern(h=%g, nu=%g) = %g", h, nu, corr);
    plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
    mxSetPr(plhs[0], &corr);
         return;
}

If I translate the C mex code above to C++ code file, I can get the C++ code compiled with the g++ complier on my mac successfully.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <boost/math/special_functions.hpp>

double matern(const double h, const double nu)
{

  double tau = sqrt(2.0*nu) * h;

  double c = 0.0;
  if(h>0.0){
    c = (pow(tau, nu) * pow(2.0, 1.0-nu) / tgamma(nu)) * boost::math::cyl_bessel_k(nu, tau);
  }else{
    c = 1.0;
  }


  return c;
}

int main(){
    double nu=0.5, h=1.0;
    double corr = matern(h, nu);
    printf("corr=%lf\n", corr);

    return 0;
}

To emphasize my question again, I do not need the C++ code, instead I would like to make my C mex code run successfully.

Answer 1

This is not really Mex-specific.

You can use the Boost special math functions from a C program; they are designed to be used in either C or C++.

To do so, you need to add

#include <boost/math/tr1.hpp>

In C++, the declarations in this file are inserted into the namespace booth::math but obviously you cannot use C++ namespaces in C code. The functions themselves are compiled with C linkage so they are global names, and can simply be used as is.

Note that Boost provides both C99 and TR1 special math functions (so it can be used with a standard C library which predates C99, if you are still living in the past), but they are found in different libraries, as explained in the Boost documentation. Even if you have a C99 standard library, after you include the Boost header you will need to use the Boost libraries for C99 functions. So you will probably need to add two link options to your build command:

-lboost_math_c99 -lboost_math_tr1

If you need float or long double versions, you will need other libraries. See the above link to the Boost documentation for details.