How does one return allocated memory from a Fortran Mex function?

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

Answer 1

[James Tursa]

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.

[Patrick]

James, thanks a ton! I had tried many variations of initializing an empty plhs(1) and using mxCalloc, but your version with mxMalloc does the trick. For the record, your "simpler approach" is obviously preferable if you can get away with it, but in my real code, I can't allocate the space until I'm already in the computational routine, so I have to do it the "hard way". Thanks for going through the trouble of the "hard way" for me!

[yuehua Feng]

Hi Patrick, I try to run your code, but not successful. The following is my code and error. My Matlab version is R2016b and Fortran version is GFORTRAN (GCC 6.3.0). Could you please help me at your convenience. Thanks very much!

#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 
mwSize :: num_i4, zero = 0
integer(4) :: mxREAL = 0
plhs(1) = mxCreateDoubleMatrix(zero,zero,mxREAL)
!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)
call mxSetM(plhs(1),num_i4)
call mxSetN(plhs(1),num_i4)
!plhs(1) = mxGetPr(something)
return
end
subroutine make_matrix(num, something, num_i4)
implicit none
mwPointer :: mxMalloc
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)
something = mxMalloc( 8 * num_i4 * num_i4 )
call mxCopyReal8ToPtr(something_r8, something, size(something_r8))
deallocate(something_r8)
return
end

The Matlab error is :

    >> mex -largeArrayDims make_matrix.F90
      Building with 'gfortran'.
      /Users/yuehuafeng/Documents/MATLAB/test/make_matrix.F90:23:48:
       call make_matrix(%val(Pr(1)), something, num_i4)
                                                      1
      Warning: Type mismatch in argument 鈥榥um鈥 at (1); passed INTEGER(8) to REAL(8)

      Error using mex
      ld: warning: object file
      (/var/folders/yv/_khzm5qs2xvf65r74s0f66j80000gn/T/mex_31931934922927_13591/make_matrix.o)
      was built for newer OSX version (10.12) than being
      linked (10.11)
      ld: warning: object file
      (/var/folders/yv/_khzm5qs2xvf65r74s0f66j80000gn/T/mex_31931934922927_13591/fortran_mexapi_version.o)
      was built for newer OSX version (10.12) than being
      linked (10.11)
      ld: warning: ignoring file
      /usr/local/lib/libgfortran.dylib, file was built for
      i386 which is not the architecture being linked
      (x86_64): /usr/local/lib/libgfortran.dylib
      Undefined symbols for architecture x86_64:
        "__gfortran_os_error", referenced from:
            _make_matrix_ in make_matrix.o
        "__gfortran_runtime_error", referenced from:
            _make_matrix_ in make_matrix.o
        "__gfortran_runtime_error_at", referenced from:
            _make_matrix_ in make_matrix.o
        "__gfortran_size0", referenced from:
            _make_matrix_ in make_matrix.o
      ld: symbol(s) not found for architecture x86_64
      clang: error: linker command failed with exit code 1
      (use -v to see invocation)

For examples in Matlab is OK:

>> mex fulltosparse.F loadsparse.F
Building with 'gfortran'.
ld: warning: object file (/var/folders/yv/_khzm5qs2xvf65r74s0f66j80000gn/T/mex_34048500661351_13591/fulltosparse.o) was built for newer OSX version (10.12) than being linked (10.11)ld: warning: 
object file (/var/folders/yv/_khzm5qs2xvf65r74s0f66j80000gn/T/mex_34048500661351_13591/loadsparse.o) was built for newer OSX version (10.12) than being linked (10.11)
ld: warning: object file (/var/folders/yv/_khzm5qs2xvf65r74s0f66j80000gn/T/mex_34048500661351_13591/fortran_mexapi_version.o) was built for newer OSX version (10.12) than being linked (10.11)
ld: warning: ignoring file /usr/local/lib/libgfortran.dylib, file was built for i386 which is not the architecture being linked (x86_64): /usr/local/lib/libgfortran.dylib

MEX completed successfully.