Fortran MEX returns only some elements of output array.
5 views (last 30 days)
Show older comments
Matthew Courtemanche
on 9 Feb 2021
Commented: Matthew Courtemanche
on 9 Feb 2021
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:
- 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?
- 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.
0 Comments
Accepted Answer
James Tursa
on 9 Feb 2021
Edited: James Tursa
on 9 Feb 2021
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.
4 Comments
James Tursa
on 9 Feb 2021
Edited: James Tursa
on 9 Feb 2021
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.
More Answers (0)
See Also
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!