MATLAB: How does one return allocated memory from a Fortran Mex function

allocatefortranMATLABmemorymex

I've been struggling to find a way to allocate an array in my Fortran Mex function and return it to Matlab. Every variation I attempt crashes, usually at the return of the Mex function. I have made a trivial example that reproduces the crux of the issue.
My first thought was just to use allocatable arrays, but that crashes upon calling the make_matrix routine. Then I tried mxCalloc, but I couldn't figure out how to make that not crash upon returning from the Mex function. My last attempts (illustrated below) were to try to use mxGetPr or mxSetPr, but any variation on that just crashes as well. I can't find an example of any of these ideas in Fortran anywhere. How does one accomplish this? It seems like it should be simple.
#include<fintrf.h>
subroutine MEXFUNCTION(nlhs, plhs, nrhs, prhs)
implicit none
MWPOINTER :: plhs(*), prhs(*)
MWPOINTER :: mxCreateDoubleMatrix, mxGetPr
MWPOINTER :: Pr(1), something
integer(kind = 4) :: nrhs, nlhs, nargin, nargout, jarg, num_i4
!Check for proper number of arguments
nargin = 1
nargout = 1
if(nrhs .ne. nargin) then
    call mexErrMsgTxt('Wrong # inputs found - expecting 1')
elseif(nlhs .ne. nargout) then
    call mexErrMsgTxt('Wrong # outputs found - expecting 1')
end if
do jarg = 1, nargin
    Pr(jarg) = mxGetPr(prhs(jarg))
end do
! Call computational subroutine:
call make_matrix(%val(Pr(1)), something, num_i4)
! Create the output array:
call mxSetPr(plhs(1), something)
!plhs(1) = mxGetPr(something)
return
end
subroutine make_matrix(num, something, num_i4)
implicit none
mwPointer :: mxCreateDoubleMatrix
real(kind = 8), intent(in) :: num
integer(kind = 4), intent(out) :: num_i4
integer(kind = 4) :: i, j
mwPointer :: something
real(kind = 8), allocatable :: something_r8(:, :)
num_i4 = int(num)
allocate(something_r8(num_i4, num_i4))
do i = 1, num_i4
    do j = 1, num_i4
        something_r8(i, j) = i * j
    end do
end do
something = mxCreateDoubleMatrix(num_i4, num_i4, 0)
call mxCopyReal8ToPtr(something_r8, something, size(something_r8))
deallocate(something_r8)
return
end

Best Answer

Fundamental Principle: You cannot mix native Fortran (or C or C++) memory within an mxArray. To do so will always result in a crash. You can only use MATLAB API function allocated memory within an mxArray. And even then, you must know what you are doing to make sure the memory is accounted for properly in the MATLAB Memory Manager so it doesn't inadvertantly free memory that is part of an mxArray, which would result in a crash.
To fix your immediate problem, note that you have used a pointer to an mxArray (the result of your mxCreateDoubleMatrix call) as a pointer to a real(8) data block. That results in writing over the mxArray structure itself in the very next line with the mxCopyReal8ToPtr call, leading to the crash. To fix that, change this line:
something = mxCreateDoubleMatrix(num_i4, num_i4, 0)
to this:
something = mxMalloc( 8 * num_i4 * num_i4 )
At that point in the code, the "something" address will be on the garbage collection list (i.e., scheduled to be free'd as soon as the mex function returns control to the caller).
Your attempt at using "something" will also cause a crash, however:
call mxSetPr(plhs(1), something)
You haven't created plhs(1) yet, so it is pointing to garbage at this point. Using it in an API call will result in a crash. To do what you are trying to do, you would need to do this:
mwSize :: num_i4, zero = 0
integer(4) :: mxREAL = 0
plhs(1) = mxCreateDoubleMatrix(zero,zero,mxREAL)
:
call mxSetPr(plhs(1),something) ! removes something from garbage collection list
call mxSetM(plhs(1),num_i4)
call mxSetN(plhs(1),num_i4)
All of the above is the hard way of doing things, however. A simpler approach is to simply use mxCreateDoubleMatrix in the calling routine to the desired size to begin with, then pass its data pointer down to the computational routine. E.g.,
mwSize :: num_i4
integer(4) :: mxREAL = 0
:
plhs(1) = mxCreateDoubleMatrix(num_i4,num_i4,mxREAL)
something = mxGetPr(plhs(1))
:
call make_matrix(%val(Pr(1)), %val(something), num_i4)
:
real*8 :: something(num_i4,num_i4)
Side Comment: You should really make the check "nlhs .ne. nargout" as "nlhs .gt. nargout" instead. The way you have it, you are not allowing the user to return the result into ans, but are forcing him/her to assign the result into a variable. That is generally too restrictive and is not how users expect most functions to operate.
Related Question