Simple newbie C mex question..

2 views (last 30 days)
Arwel
Arwel on 12 Jul 2016
Edited: James Tursa on 8 Aug 2016
Hi, I'm (very) new to C and need to write a C++ mex file that inputs and outputs some cell arrays. So, I'm trying to get my head around what I need to do to transfer and manipulate them in the correct format and so on. So, below is my simple 'noddy' test program, that will accept a cell array and then tell me how many dimensions it has. I then try to pass the same array to a small function to do the same thing. This compiles, but when I run it I get the following....
>> cc = {[1,2,3],[4,5,6],[7 8 9]};
>> gateWayTest2(cc)
Inside gateway , var has 3 dims
Inside function, var has 1 dims
>>
So, I would expect 3 dims in both. What am I doing wrong here? (I can partially answer my own question in that I'm very new to C and am struggling to get my head around pointers...). But anyway, the stage I want to get to is that I can do some maths on the contents of 'cc' in the function, and then output this modified 'cc' back into Matlab. So if anyone could show me how I could do something simple like add a constant to some of the elements and then output that, that would be a really really useful! Many thanks in advance, Arwel
#include "mex.h"
void testFunction(mxArray *in)
{
double *numberOfContrasts;
numberOfContrasts = mxGetNumberOfElements(in);
mexPrintf("Inside function, var has %d dims \n",numberOfContrasts);
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
double *numberOfContrasts;
mxArray *xdata;
/*Check for correct number and type of input arguments*/
if (nrhs != 1) {
mexErrMsgIdAndTxt("MATLAB:cellTest:invalidNumInputs",
"One input argument required");
};
if (mxIsCel!=1) {
mexErrMsgIdAndTxt("MATLAB:cellTest:invalidType",
"Cell array required");
};
/*Print how many dimensions in cell array as a test*/
numberOfContrasts = mxGetNumberOfElements(prhs[0]);
mexPrintf("Inside gateway , var has %d dims \n",numberOfContrasts);
/*Pass to a function to a simple fn make sure it 'arrives' in the correct format*/
xdata = mxGetPr(prhs[0]);
testFunction(xdata);
return;
}

Answers (5)

James Tursa
James Tursa on 12 Jul 2016
Edited: James Tursa on 12 Jul 2016
A few problems:
1) This code is wrong:
double *numberOfContrasts;
numberOfContrasts = mxGetNumberOfElements(in);
mexPrintf("Inside function, var has %d dims \n",numberOfContrasts);
The return value of mxGetNumberOfElements is size_t, but you are assigning that return value into a pointer_to_double. Then you print the pointer with a %d format which is intended for integer types. There is no need to use a pointer here. Simply use the same type as the return value. E.g.,
size_t numberOfContrasts; // <-- changed type to size_t
numberOfContrasts = mxGetNumberOfElements(in);
mexPrintf("Inside function, var has %d dims \n",numberOfContrasts);
2) Same problem in your mexFunction code:
double *numberOfContrasts; // <-- Need to change (double *) to size_t
:
numberOfContrasts = mxGetNumberOfElements(prhs[0]);
mexPrintf("Inside gateway , var has %d dims \n",numberOfContrasts);
3) You are passing a data pointer to a routine that expects an mxArray pointer:
xdata = mxGetPr(prhs[0]);
testFunction(xdata);
In the above code, the data pointer for prhs[0] points to a memory block that contains three mxArray addresses. It does not point to memory block that contains an mxArray. (In fact, I would not have been surprised if this code had crashed MATLAB with an invalid memory access violation.) The correct way to retrieve the address of this data block would be:
mxArray **cellblock;
:
cellblock = (mxArray **) mxGetData(prhs[0]);
That is, cellblock is a pointer, and it points to memory that contain (mxArray *) types.
All that being said, you don't really want to be doing this anyway in your code. All you want is to pass prhs[0] directly. E.g.,
testFunction(prhs[0]);
This should generate a compiler warning since you are passing a const value to a routine that is not expecting a const value. You can correct that by doing this:
void testFunction(const mxArray *in)
4) I assume this is a typo:
if (mxIsCel!=1) {
and you really have this in your code:
if (mxIsCell(prhs[0]) != 1) {
5) Finally, I would strongly advise that when you mex your code you use the verbose option -v. That way a lot of these programming errors that may only appear as warnings to the compiler will be shown on the screen to you.

Arwel
Arwel on 12 Jul 2016
Thanks James for your help.... I am a Matlab coder, not a C one, so I'm really not used to thinking about pointers. I've been thinking about this on the way home tonight and I'm a bit clearer about what I would really like to do. This is... (1) My data in matlab is a simple cell array of doubles. This I want to pass to C.. (2) In the wrapper function, I want to de-reference all this right down to something that resembles a cell array that I would recognise as a Matlab coder. (3) Pass this to the calculation function. (4) In my calculation function work directly on this array in a 'matlaby' way (i.e. no pointers). The return is another cell array of doubles. (5) clean up, and pack this cell of doubles in the correct way to send back to Matlab.
Is this even possible or am I hoping against hope in this pointers of pointers or whatever wasteland? If I'm honest I don't really care whats actually going on in terms what is a pointer and so on, just so long as when it's in my function I can just do the maths on the elements of this and just get on with things. If this is possible do you think you could make me a very simple 'hello world'? Many Thakns, Arwel
  3 Comments
Arwel
Arwel on 13 Jul 2016
Usual answer - Speed! I'm not a complete newbie to mex files but my background in this is with Fortran mex, which for whatever reason make far more sense to me. But, I want to go even faster and I want to use OpenMP within the mex on some loops (yes... I KNOW about the parallel toolbox, but this is going to be compiled code which is going to be shared gratis to academics, and deployed so they don't have to buy Matlab - so the 'you need the Distributed Toolbox....' blah blah that the Parallel toolbox requires is totally and completely useless for me). I need threads not 'workers'.... Anyway, I know that OpenMP works with Fortran too, but the limit then is the really small number of compilers that will work with Matlab. If you google around there are plenty of examples of using OpenMP in C++ mex files (even with GCC..), but not a single one with Fortran. So I'll have to learn to do this with C....
James Tursa
James Tursa on 13 Jul 2016
OK. So can you elaborate just a bit on the calculations you will be doing? Are you thinking of manipulating a cell array input variable in-place? And multi-threading that work out? Both of these are not easy in a mex routine for two main reasons:
(1) It can be very difficult to tell when a cell element is shared with another variable in a mex routine, but you need to know this in order to manipulate the variable properly. There is really nothing in the API that you can use for this ... you have to hack into the mxArray variable structure. Not pretty.
(2) None of the API function that allocate memory are thread-safe. So if you were thinking of multi-threading anything that involved creating cell elements you can forget that right now ... it would likely mess up the MATLAB Memory Manager and crash MATLAB. The only practical solution is to create all of the cell elements up front, before you enter the parallel section, and then in the parallel section you can fill in the data.

Sign in to comment.


Arwel
Arwel on 19 Jul 2016
Hi James, Well, I switched back to Fortran and I thought I had it worked out, but unfortunately I'm still not getting anywhere. At the moment my fortran MEX is called from within a loop in matlab, with the inputs either being double arrays or scalars, like so....
C=======================================================================
C Abeles Calculation for Matlab
C
C out = abeles(xdata,sld,nbair,nbsub,ssub,repeats,s_between_repeats,resol)
C
C=======================================================================
#include "fintrf.h"
subroutine mexfunction(nlhs, plhs, nrhs, prhs)
mwPointer, external :: mxGetPr, mxCreateDoubleMatrix
mwSize, external :: mxGetM, mxGetN
integer nlhs, nrhs
mwPointer plhs(*), prhs(*)
mwPointer pr_in, pr_in2, pr_in3, pr_in4, pr_in5, pr_in6, pr_in7
mwPointer pr_out, pr_in8
integer m_in_q, n_in_q
integer size, points, layers
integer m_in_sld, n_in_sld
if(nrhs .ne. 8) then
call mexErrMsgTxt('8 inputs required to Abeles.')
endif
if(nlhs .gt. 1) then
call mexErrMsgTxt('Less than one output required to Abeles.')
endif
! work out points
m_in_q = mxGetM(prhs(1))
n_in_q = mxGetN(prhs(1))
if (n_in_q .gt. 1) then
call mexErrMsgTxt('Qz needs to be a column vector')
endif
points = m_in_q
! do the same for layers
m_in_sld = mxGetM(prhs(2))
n_in_sld = mxGetN(prhs(2))
if (n_in_sld .ne. 3) then
call mexErrMsgTxt('SLD needs three columns: d r s')
endif
layers = m_in_sld
pr_in = mxGetPr(prhs(1))
pr_in2 = mxGetPr(prhs(2))
pr_in3 = mxGetPr(prhs(3))
pr_in4 = mxGetPr(prhs(4))
pr_in5 = mxGetPr(prhs(5))
pr_in6 = mxGetPr(prhs(6))
pr_in7 = mxGetPr(prhs(7))
pr_in8 = mxGetPr(prhs(8))
plhs(1) = mxCreateDoubleMatrix(m_in_q, n_in_q, 0)
pr_out = mxGetPr(plhs(1))
C Call the computational routine.
call comp(%val(pr_out),%val(pr_in),%val(pr_in2),%val(pr_in3),
+ %val(pr_in4),%val(pr_in5),%val(pr_in6),
+ %val(pr_in7),%val(pr_in8),layers,points)
return
end
C=========================================================================
c
C computational routine.....
C
C=========================================================================
subroutine comp(out,x,sld,nbair,nbsub,ssub,nrepeats,resType,
+ resol,layers,points)
integer size, i
integer layers,points,j,reploop
double complex c0 , c1 , ci , beta , rij , N , M , MI , rimaj
double complex blast , num , den , quo
double complex pimag , pair , psub , plast
real*8 pi , twopi , loop , nl , lam , nrepeats, ilow, ihi
real*8 thick , q , theta, rho, rough, resol
real*8 preal , snair , snsub , snlay, rfac, resType
real*8 out(points,1),x(points,1),sld(layers,3),nbair,nbsub,ssub
real*8 dummydata(points,1), g(points,1)
dimension N(2,2) , M(2,2) , MI(2,2)
..etc
In this 'xdata' and 'sld' are double arrays, and the rest are scalar doubles. This is called from within a loop in Matlab. What I want to be able to do is call this once from matlab, passing it all the data at once. So it would get xdata as {xdata1 ; xdata2 ; xdata3 ..... xdata_n}, the same for SLD, the others now as arrays, and then just call comp in a loop. Then, the output would be a single cell array of the n results of the call to comp.
However, I'm getting absolutely nowhere with this and just can't make it work :( Could you please, please show me what to do here, because this is driving me nuts now.... Cheers, Arwel
  3 Comments
Arwel
Arwel on 4 Aug 2016
Hi James, I've been distracted by some other things, but have just come back to this today. I've got this far with my Fortran....
#include "fintrf.h"
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use omp_lib
C Declarations
implicit none
C Mex function arguments
mwPointer plhs(*), prhs(*), dat
integer nlhs, nrhs
C Function Declarations
mwPointer mxCreateCellMatrix
mwPointer mxGetCell
mwPointer mxGetData
mwPointer mxGetPr
mwPointer mxGetNumberOfElements
mwPointer mxCreateDoubleMatrix
mwPointer mxGetM, mxGetN
mwPointer mxDuplicatearray
mwSize i, m, n
integer*4 :: ComplexFlag = 0
integer*4 :: layers, points
real*8 :: test(1,1)
C Pointers Arrays and vars
mwPointer numberOfContrasts
mwPointer outputArray
mwPointer inputArray_xdat, inputArray_sld, inputArray_nbair
mwPointer inputArray_nbsub, inputArray_ssub, inputArray_rep
mwPointer inputarray_resol
mwPointer thisCell_xdat, thisCell_sld
mwPointer calcOutpArray
mwPointer thisArray, thisArray2
mwPointer size
mwPointer pr_out, pr_in_xdata
mwPointer pr_in_sld, pr_in_nbair, pr_in_nbsub
mwPointer pr_in_rep, pr_in_resol, pr_in_ssub
C input checking...
if(nrhs .lt. 7) then
call mexErrMsgTxt('7 inputs required to Abeles.')
endif
if(nlhs .gt. 1) then
call mexErrMsgTxt('One output required from Abeles.')
endif
C Get the number of conts from size of the xdata cell array..
inputArray_xdat = prhs(1)
numberOfContrasts = mxGetNumberOfElements(inputArray_xdat)
inputArray_sld = prhs(2)
inputArray_nbair = prhs(3)
inputArray_nbsub = prhs(4)
inputArray_rep = prhs(5)
inputArray_resol = prhs(6)
inputArray_ssub = prhs(7)
pr_in_nbair = mxGetPr(inputArray_nbair)
pr_in_nbsub = mxGetPr(inputArray_nbsub)
pr_in_rep = mxGetPr(inputArray_rep)
pr_in_resol = mxGetPr(inputArray_resol)
pr_in_ssub = mxGetPr(inputArray_ssub)
C Output cell array will have the same number of contrasts
m = numberOfContrasts
n = 1
plhs(1) = mxCreateCellMatrix(m,n)
outputArray = mxCreateCellMatrix(m,n)
c Loop over all the elements in the input array and call comp...
!$OMP PARALLEL
!$OMP DO
do i=1,numberOfContrasts
c Pointers for input variables
thisCell_xdat = mxGetCell(inputArray_xdat,i)
points = mxGetN(thisCell_xdat)
pr_in_xdata = mxGetPr(thisCell_xdat)
thisCell_sld = mxGetCell(inputArray_sld,i)
layers = mxGetM(thisCell_sld)
pr_in_sld = mxGetPr(thisCell_sld)
call mxCopyPtrToReal8(pr_in_nbair,test,3)
c Make the array for output from calcfn
n = 1
calcOutpArray = mxCreateDoubleMatrix(points,n,ComplexFlag)
pr_out = mxGetPr(calcOutpArray)
call abeles(%VAL(pr_out),%VAL(pr_in_xdata),%VAL(pr_in_sld),
+ %VAL(pr_in_nbair),%VAL(pr_in_nbsub),%VAL(pr_in_rep),
+ %VAL(pr_in_resol),%VAL(pr_in_ssub),
+ layers,points,i,numberOfContrasts)
call mxSetCell(outputArray,i,calcOutpArray)
10 end do
!$OMP END DO
!$OMP END PARALLEL
plhs(1) = mxDuplicatearray(outputArray)
return
end
Now, if I compile this with
>> mex -v abeles_loop.f
.. it works just fine. But if I go...
>> setenv('OMP_NUM_THREADS','4') >> mex -v -g COMPFLAGS='$COMPFLAGS -openmp' abeles_loop.f
..then it seg faults :(
What am I missing here?? I'm worried about using the mxSetCell inside the parallel loop like that maybe?
James Tursa
James Tursa on 4 Aug 2016
You can't do anything inside of a parallel loop that uses the MATLAB Memory Manager to allocate memory. So this line is trouble:
calcOutpArray = mxCreateDoubleMatrix(points,n,ComplexFlag)
The other API calls inside your parallel loop (mxGetCell, mxGetM, mxGetN, mxGetPr) do not allocate memory so they should be OK in a parallel loop.
The mxSetCell call should also be OK since I don't think it will allocate memory, but I am not 100% sure about this. To my knowledge, this line:
call mxSetCell(outputArray,i,calcOutpArray)
does the following:
- Puts the calcOutpArray mxArray address into the i'th spot of outputArray.
- Changes the type of calcOutpArray to SUB_ELEMENT (this is done in-place, no memory allocation required)
- Removes the calcOutpArray mxArray address from the garbage collection list (again, I don't see how this act could result in memory allocation)
But the reason I am not 100% sure about memory allocation is that since R2011b MATLAB has included something in the pi data area of cell arrays. That pointer used to be NULL for cell arrays, but now it has something in it and I have no idea what that is. If there is something in the background there that allocates memory for some reason, then mxSetCell would be a problem in a parallel loop. (If that turns out to be the case, it is an easy workaround to just save the calcOutpArray mxArray pointers off to the side inside the parallel loop and then use mxSetCell on all of them after the parallel loop)
BOTTOM LINE:
You need to rewrite this parallel loop so that there is no MATLAB Memory Manager allocations inside of it. You basically need to do all of the anticipated mxCreateDoubleMatrix etc calls ahead of time and then just fill in the data inside the parallel loop.

Sign in to comment.


Arwel
Arwel on 5 Aug 2016
(This is with R2015b on Windows 7, with Visual Studio 2013 by the way)

Arwel
Arwel on 5 Aug 2016
Hi James, Okay..... I got something that actually runs with my test problem (6.33 seconds single threaded, 2.12 seconds OMP! :-), but it's not quite there still. Plugging it into the real calculations it looks initially like its working, but then soon starts to misbehave with some strange effects appearing. Taking the mxCreateDoubleMatrix out of the loop definitely worked, but it's still not quite there. Here's what I tried...
inputArray_xdat = prhs(1)
numberOfContrasts = mxGetNumberOfElements(inputArray_xdat)
inputArray_sld = prhs(2)
inputArray_nbair = prhs(3)
inputArray_nbsub = prhs(4)
inputArray_rep = prhs(5)
inputArray_resol = prhs(6)
inputArray_ssub = prhs(7)
pr_in_nbair = mxGetPr(inputArray_nbair)
pr_in_nbsub = mxGetPr(inputArray_nbsub)
pr_in_rep = mxGetPr(inputArray_rep)
pr_in_resol = mxGetPr(inputArray_resol)
pr_in_ssub = mxGetPr(inputArray_ssub)
C Output cell array will have the same number of contrasts
m = numberOfContrasts
n = 1
plhs(1) = mxCreateCellMatrix(m,n)
outputArray = mxCreateCellMatrix(m,n)
c Fill the output array with double arrays of correct size
n = 1
do i=1,numberOfContrasts
thisCell_xdat = mxGetCell(inputArray_xdat,i)
points = mxGetN(thisCell_xdat)
calcOutpArray = mxCreateDoubleMatrix(points,n,ComplexFlag)
call mxSetCell(outputArray,i,calcOutpArray)
end do
c Loop over all the elements in the input array and call comp...
!$OMP PARALLEL
!$OMP DO
do i=1,numberOfContrasts
c Pointers for input variables
thisCell_xdat = mxGetCell(inputArray_xdat,i)
points = mxGetN(thisCell_xdat)
pr_in_xdata = mxGetPr(thisCell_xdat)
thisCell_sld = mxGetCell(inputArray_sld,i)
layers = mxGetM(thisCell_sld)
pr_in_sld = mxGetPr(thisCell_sld)
c Get the array for output from calcfn
calcOutpArray = mxGetCell(outputArray,i)
pr_out = mxGetPr(calcOutpArray)
call abeles(%VAL(pr_out),%VAL(pr_in_xdata),%VAL(pr_in_sld),
+ %VAL(pr_in_nbair),%VAL(pr_in_nbsub),%VAL(pr_in_rep),
+ %VAL(pr_in_resol),%VAL(pr_in_ssub),
+ layers,points,i,numberOfContrasts)
10 end do
!$OMP END DO
!$OMP END PARALLEL
plhs(1) = mxDuplicatearray(outputArray)
return
end
  1 Comment
James Tursa
James Tursa on 8 Aug 2016
Edited: James Tursa on 8 Aug 2016
I think you need to make your scalar variables private in the parallel loop. By default, I don't think OpenMP makes these variables private. So try adding a private clause to your OpenMP meta command and make the following variables private: thisCell_Xdat, points, pr_in_xdata, thisCell_sld, layers, pr_in_sld, calcOutpArray, pr_out. I think the iteration index variable i is the only variable that is private by default.
Also, these lines are doing extra work that is not needed (and in fact create a temporary memory leak):
plhs(1) = mxCreateCellMatrix(m,n)
outputArray = mxCreateCellMatrix(m,n)
:
plhs(1) = mxDuplicatearray(outputArray)
You have created three separate cell matrices, two of which get garbage collected (otherwise they would have leaked). You only need one cell matrix here. E.g.,
plhs(1) = mxCreateCellMatrix(m,n)
outputArray = plhs(1)
:
! plhs(1) = mxDuplicatearray(outputArray) <-- Get rid of this line

Sign in to comment.

Categories

Find more on Fortran with MATLAB in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!