MATLAB: Fortran MEX returns only some elements of output array.

fortranmex

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.

Best Answer

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.
Related Question