I doubt you will be able to do this. Inside of mex routines, the API functions that deal with memory management are not (or at least they have not been in the past) thread safe. E.g., all of the mxCreateETC functions dynamically allocate memory for a new mxArray, so those functions are not thread safe. And the memory allocation functions mxMalloc, mxCalloc, and mxRealloc are not thread safe because they allocate memory. mxGetProperty also allocates memory and is not thread safe. Several other functions are in this category as well. So unless things have fundamentally changed in the API function libraries recently, you can't do anything that involves the MATLAB Memory Manager in parallel threads safely.
Functions that are thread safe would be functions like mxGetPr, mxGetPi, mxIsDouble, etc ... i.e., functions that simply return a value or status and do not allocate memory. Those could be safely used inside parallel threads.
Unfortunately, mexCallMATLAB allocates memory for the return value in the plhs array, so mexCallMATLAB is not thread safe in general. Even if mexCallMATLAB did not return anything (e.g., nlhs = 0), it would not surprise me that the memory allocations going on in the MATLAB function you called (i.e., CALC) were not thread safe.
My best advice would be to convert CALC into C code and then put that inside your OpenMP threads. If there is memory allocation involved, do it outside the parallel code.
You could of course ignore all of the above and simply code your CALC call inside the loop and see what happens, but I suspect it will crash MATLAB at some point. Consider this line:
D(:, i) = CALC(A, B(:,i), C(i));
The bulk of the C code for the above line would be getting the mxArray's for B(:,i) and C(i), and then copying the result from mexCallMATLAB into D(:,i). An outline of of how you would do this is as follows:
#include <string.h>
size_t m;
double *d, *b;
mxArray *Ci, *Bi;
mxArray rhs[3], lhs[1];
:
m = mxGetM(prhs[1]);
:
for( i=0; i<n; i++ ) {
:
Ci = mxCreateDoubleScalar( *C++ );
Bi = mxCreateDoubleMatrix( m, 1, mxREAL );
b = mxGetPr(Bi);
memcpy( b, B, m * sizeof(*b) );
B += m;
rhs[0] = prhs[0];
rhs[1] = Bi;
rhs[2] = Ci;
mexCallMATLAB( 1, lhs, 3, rhs, "CALC" );
mxDestroyArray(Bi);
mxDestroyArray(Ci);
d = mxGetPr(lhs[0]);
memcpy( D, d, r * sizeof(*d) );
D += r;
mxDestroyArray(lhs[0]);
}
