Fortran MEX returns only some elements of output array.

Hi all, I'm trying to learn how to write MEX functions using Fortran. I have encountered a couple issues, and the documentation is quite dated (mostly fixed-form Fortran, and predating the current -R2018a API). Many suggestions lead to functions like mxGetPr and mxGetPi which are 'not recommended' in the current API. This goes for all the .F examples in $matlabroot/extern/examples/.

I am on Ubuntu 20.04, Matlab R2019b, gfortran 9.3.0

I have two issues:

1. If I declare my arrays as allocatable, and then allocate them, I will hit a segfault or some other instability causing Matlab to lock or crash. The "solution" is to pre-allocate the arrays with some maximum size, and then work with inputs smaller than that hidden limit. It no longer segfaults, but this harms the ability to write a general-purpose array input without adding an arbitrarily large allocation for all arrays. I'm not sure of a way to use `%VAL(A)` with matrix operations -- perhaps `RESHAPE()` into a local variable that is matrix-shaped within the computational subroutine?
2. In my minimum-working training example, I just pass array A to a subroutine, and then copy it to B, and return B. Some, but not all of the elements of A are returned in B. The usage in Matlab is B = copyab(A).

The F90 file, 'copyab.F90' (both gateway and computational subroutines)

#include "fintrf.h"
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
  implicit none

  mwPointer :: plhs(*), prhs(*)
  integer :: nlhs, nrhs

  mwPointer :: mxGetDoubles
  mwPointer :: mxCreateDoubleMatrix
  mwPointer :: mxGetM, mxGetN
  mwPointer :: A_ptr, B_ptr
  mwPointer :: m, n
  mwSize :: mn

  ! Issue here: using allocatable, allocate() causes segfault.
  integer, parameter :: maxm = 10 
  integer, parameter :: maxn = 10
  real :: A(maxm,maxn), B(maxm,maxn)
!  real, allocatable :: A(:,:), B(:,:)
  m = mxGetM(prhs(1))
  n = mxGetN(prhs(1))
  mn = m * n
!  allocate(A(m,n), B(m,n))
  
  ! RHS: input args
  A_ptr = mxGetDoubles(prhs(1))
  call mxCopyPtrToReal8(A_ptr, A, mn)
  
  ! LHS: return vals
  plhs(1) = mxCreateDoubleMatrix(m, n, 0)
  B_ptr = mxGetDoubles(plhs(1))
  call copyab(A, B, m, n)
  call mxCopyReal8ToPtr(B, B_ptr, mn)
end subroutine mexFunction


subroutine copyab(A, B, m, n)
  implicit none

  mwSize :: m, n
  real :: A(m,n), B(m,n)
  
  B = A

end subroutine copyab

Makefile with compile options:

FC = gfortran
FCFLAGS = -O3 -cpp -fPIC -ffree-form -fdefault-real-8 -fdefault-integer-8 -Wall
MATLABDIR = /usr/local/MATLAB/R2019b
MEX = $(MATLABDIR)/bin/mex
MEXLIBDIR = $(MATLABDIR)/sys/os/glnxa64
MEXLIBS = -lgfortran
MEXINCLUDE = $(MATLABDIR)/extern/include
MEXFLAGS = -R2018a COMPFLAGS='$(FCFLAGS)'
EXT = mexa64
copyab.$(EXT): copyab.F90
  $(MEX) $^ $(MEXFLAGS) -output $@ -I$(MEXINCLUDE) -L$(MEXLIBDIR) $(MEXLIBS)

clean:
  rm -rf *.$(EXT) *.o *.mod

And the test Matlab code, 'test.m'

A = rand(3,3)
B = copyab(A)

which results in,

>> test

A =

    0.8147    0.9134    0.2785
    0.9058    0.6324    0.5469
    0.1270    0.0975    0.9575


B =

    0.8147    0.9134    0.0000
    0.9058    0.0000    0.0000
    0.1270    0.0000    0.0000

B appears to containg floor(numel(A)/2) correct elements, regardess of the shape of A.

Answer 1

[James Tursa]

The problem isn't using allocatable variables ... MATLAB mex functions can handle that just fine. The problem is that you are using two different types of floating point variables, single (aka real) and double (aka real*8). You allocate the variable as real (4-bytes per element) but then copy it as if it is real*8 (8 bytes per element), so you access invalid memory and get the seg fault. The fix is to be consistent with your floating point types ... use either real or real*8. If you use real then you would need to use the mxCopyPtrToReal4 and mxCopyReal4ToPtr functions. If you use real*8 then you can use the mxCopyPtrToReal8 and mxCopyReal8ToPtr functions. Also for readability you should be using mxGetSingles for real variables, although this doesn't really make a difference in Fortran since the address is just returned in an integer variable. You could also be generic and just use mxGetData for this, which will work in both the R2017b and R2018a API memory models.

I would also strongly advise putting in argument checking code in mex functions to avoid unexpected arguments causing erroneous results or seg faults.

As an aside, what you are doing is fine for small variables where the deep data copies don't hurt performance much. But if you are using very large variables in your actual problem, you would want to avoid the deep data copies if at all possible. That would mean doing something like using the %VAL(etc) construct in your code to pass pointers by value to implicit interface subroutines. E.g., you could do something like this (Caveat ... I don't have a Fortran compiler so this is all untested):

#include "fintrf.h"
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
  implicit none

  mwPointer :: plhs(*), prhs(*)
  integer :: nlhs, nrhs

  ! using external makes it clear to the reader that these are functions
  mwPointer, external :: mxGetDoubles
  mwPointer, external :: mxCreateDoubleMatrix
  mwPointer, external :: mxGetM, mxGetN
  integer*4, external :: mxIsDouble, mxIsComplex, mxIsSparse
  mwPointer :: A_ptr, B_ptr
  mwSize :: m, n ! changed to mwSize to match mxCreateDoubleMatrix signature exactly
  integer*4 :: ComplexFlag = 0 ! added explicit type to match mxCreateDoubleMatrix signature

  ! adding argument checks
  if( nlhs > 1 ) then
      call mexErrMsgTxt('Too many outputs')
  endif
  if( nrhs != 1 ) then
      call mexErrMsgTxt('Need one full real double input')
  endif
  if( mxIsDouble(prhs(1))==0 .or. mxIsComplex(prhs(1))==1 .or. mxIsSparse(prhs(1))==1 ) then
      call mexErrMsgTxt('Need one full real double input')
  endif
  
  m = mxGetM(prhs(1)) ! mwSize = mwPointer
  n = mxGetN(prhs(1)) ! mwSize = mwPointer
  
  ! RHS: input args
  A_ptr = mxGetDoubles(prhs(1))
  
  ! LHS: return vals
  plhs(1) = mxCreateDoubleMatrix(m, n, ComplexFlag) ! using explicit typed variable for argument
  if( m*n /= 0 ) then  ! check for empty variable
      B_ptr = mxGetDoubles(plhs(1))
      call copyab(%VAL(A_ptr), %VAL(B_ptr), m, n) ! pass addresses by value
  endif
end subroutine mexFunction


subroutine copyab(A, B, m, n)
  implicit none

  mwSize :: m, n
  real*8 :: A(m,n), B(m,n) ! changed to real*8 to match doubles in calling routine exactly
  
  B = A

end subroutine copyab

Note that you should never pass literal integers to Fortran API functions, since in many cases you can't guarantee what the size of integer the compiler will use (depends on mex and compiler settings). Always use variables that are explicitly typed according to the API interface documentation to avoid problems.

[Matthew Courtemanche]

Thanks for the tips.

This lead me to find the issue was that passing the gfortran option -ffdefault-real-8 to COMPFLAGS in mex was NOT compiling with -ffdefault-real-8, and subsequently the reals were left as real(4)s. If I explicitly set real(8) within the computational subroutine, it works correctly. I'll need to figure out how to get the gfortran options to work.

As well, the way you have %VAL seems to work(), so I can switch over to that.

[Matthew Courtemanche]

For the record, passing

mex ... FFLAGS='-fdefault-real-8 -fdefault-integer-8'

was the correct syntax, and -fdefault-real-8, etc. is respected. COMPFLAGS is not correct.

Apparently calling mex -setup Fortran should point to a default configuration .sh file, but this is not the case for me.

[James Tursa]

Did you get rid of the silly /fixed option that TMW includes in their Fortran build files? Although I have complained about this to TMW in the past, as of R2020a it is still there (at least for the Intel Fortran xml build files) and I have no idea why.

[Matthew Courtemanche]

I did not change the .xml, I passed `-ffree-form` into FFLAGS, and that seems to have worked. I suppose this means I will need to do this manual override of the fixed form default in every Makefile.

I have not used an Intel compiler for Fortran-MEX yet (only Linux w/GCC), but I will probably need to build my program on all platforms later on, and if can do an ifort equivalent to gfortran -ffree-form I shall report back.