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.

I think you've said you can compute the result you want with C++ code that you know how to produce. So what's the problem? Surely not as simple as needing to declare such a function to have C linkage? — John Bollinger, Jun 25 '18 at 15:21
Hi @JohnBollinger, I have elaborated my post. What I would like to do is to make my C mex code run successfully. — Bayes, Jun 25 '18 at 17:38
Hi @rici, thanks for your suggestion. I tried to include that in my C mex code, but it still does not work. Any further suggestion? — Bayes, Jun 25 '18 at 17:39
I still do not understand. You appear already to have all the pieces you need. Do you in fact know what "C linkage" means in C++? Do you not see how that might be applicable to your problem? Is there some reason why it would not be applicable? — John Bollinger, Jun 25 '18 at 18:17
Hi @JohnBollinger, as far as I know, MATLAB does not support support such functionality, given the fact that there is no such official documentation from MathWork. Maybe an export with knowledge of C/C++ Mex functions in MATLAB can give a definitive answer to this. — Bayes, Jun 25 '18 at 18:50
@rici Thanks for your suggestion. After I tried what you suggested, I got the following error message in MATLAB: Error using mex Undefined symbols for architecture x86_64: "_boost_tgamma", referenced from: _matern in matern.o _mexFunction in matern.o ld: symbol(s) not found for architecture x86_64 clang: error: linker command failed with exit code 1 (use -v to see invocation) This message shows that I have trouble to use tgamma function. I thought tgamma function can be directly accessed from math.h. — Bayes, Jun 26 '18 at 13:20
@rici Thanks for that. I think I linked the boost library correctly, since 1) if I delete the tgamma function but keep the boost_cyl_bessel_k, the code will be compiled successfully. 2) if I delete boost_cyl_bessel_k and the boost include line but keep the tgamma function, the code will be compilied successfully. BTW, the compilation code I used is mex -v -L/usr/local/lib -lboost_math_tr1 -I/usr/local/include/ matern.c. — Bayes, Jun 26 '18 at 13:51
@rici Thank you very much. It works now. Just one more question on this. Why do I need to include the -lboost_math_c99 in order to use c99 fucntions in the option code after I use the boost library (i.e., TR1 functions)? I thought the math.h provides tgamma function. In my code, I use tgamma instead of boost_tgamma. I — Bayes, Jun 26 '18 at 15:05
@pulong: Because the boost header substitutes its own functions (using macros). I don't know why it does that, but it might be explained somewhere in the Boost documentation. The page I referenced does say that it does that. — rici, Jun 26 '18 at 15:07
@pulong: I cleaned up my comments and wrote an answer, since we now know what you needed and it may be of use to other people. — rici, Jun 26 '18 at 15:16

score 0 · Accepted Answer · answered Jun 26 '18 at 15:15

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.

Using Besselk function in C MexFunction

1 Answers1