MATLAB: Matlab crashes while running MEX Fortran 90 routine

mex fortran

I compiled a Fortran mex file based off someone's code in Matlab 2019a using Intel fortran compiler 2019 and Microsoft visual studio 2017. Everytime Matlab uses the mex file to run the FortraN routine, it crashes. I am not sure if there is some error in the way I compiled the mex file or some issue with the code.
To compile the mex file, I used the command "mex('-largeArrayDims','-silent','-f','intel_fortran_19_vs2017.xml', 'mvndstpack.F90')" in Matlab 2019a (I have . Here is my code:
#include "fintrf.h"
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
    ! Declarations
	implicit none
    ! mexFunction argument
    mwPointer plhs(*), prhs(*)    
    integer*4 nlhs, nrhs   
    ! Function declarations
    mwSize mxGetN  
    mwpointer mxGetPr, mxCreateNumericArray, mxGetDimensions 
    double precision mxGetScalar  
    integer*4 mxClassIDFromClassName     
    ! Pointers to input/output mxArrays
    mwpointer x8_pr,x9_pr
    mwpointer pf1_pr,pf2_pr,pf3_pr
    mwpointer o1_pr,o2_pr,o3_pr
    
    ! Array information
	mwSize nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9,nodes
	integer*4 myclassid
	double precision, allocatable, dimension(:) :: x8,x9
	double precision x8i,x9i
	double precision, allocatable, dimension(:,:,:,:,:,:,:,:,:) :: pf1,pf2,pf3
	double precision, allocatable, dimension(:,:,:,:,:,:,:) :: o1,o2,o3
	  
    ! Load Inputs
    nx1 = mxGetScalar(prhs(1))
    nx2 = mxGetScalar(prhs(2))
    nx3 = mxGetScalar(prhs(3))
    nx4 = mxGetScalar(prhs(4))
    nx5 = mxGetScalar(prhs(5))
    nx6 = mxGetScalar(prhs(6))
    nx7 = mxGetScalar(prhs(7))
    nx8 = mxGetScalar(prhs(8))
    nx9 = mxGetScalar(prhs(9))
    nodes = nx1*nx2*nx3*nx4*nx5*nx6*nx7*nx8*nx9
    allocate(x8(nx8))
    allocate(x9(nx9))
    x8_pr = mxGetPr(prhs(10))
    x9_pr = mxGetPr(prhs(11))
    call mxCopyPtrToReal8(x8_pr,x8,nx8)
    call mxCopyPtrToReal8(x9_pr,x9,nx9)
    ! Point to evaluate
    x8i = mxGetScalar(prhs(12))
    x9i = mxGetScalar(prhs(13))
	  
    ! Policy functions
    allocate(pf1(nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9))
    allocate(pf2(nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9))
    allocate(pf3(nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9))
    pf1_pr = mxGetPr(prhs(14))
    pf2_pr = mxGetPr(prhs(15))
    pf3_pr = mxGetPr(prhs(16))
    call mxCopyPtrToReal8(pf1_pr,pf1,nodes)
    call mxCopyPtrToReal8(pf2_pr,pf2,nodes)
    call mxCopyPtrToReal8(pf3_pr,pf3,nodes)
	  
    !Create array for return argument
    myclassid = mxClassIDFromClassName('double')
    allocate(o1(nx1,nx2,nx3,nx4,nx5,nx6,nx7))
    allocate(o2(nx1,nx2,nx3,nx4,nx5,nx6,nx7))
    allocate(o3(nx1,nx2,nx3,nx4,nx5,nx6,nx7))
    plhs(1) = mxCreateNumericArray(7,[nx1,nx2,nx3,nx4,nx5,nx6,nx7],myclassid,0)
    plhs(2) = mxCreateNumericArray(7,[nx1,nx2,nx3,nx4,nx5,nx6,nx7],myclassid,0)
    plhs(3) = mxCreateNumericArray(7,[nx1,nx2,nx3,nx4,nx5,nx6,nx7],myclassid,0)
    o1_pr = mxGetPr(plhs(1))
    o2_pr = mxGetPr(plhs(2))
    o3_pr = mxGetPr(plhs(3))
                      
    ! Call subroutine for assignment
    call allterp(	nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9, &
                                        x8,x9, &
					x8i,x9i, &
					pf1,pf2,pf3, &
					o1,o2,o3)
    
    ! Load Fortran array to pointer (output to MATLAB)
    call mxCopyReal8ToPtr(o1,o1_pr,nx1*nx2*nx3*nx4*nx5*nx6*nx7)
    call mxCopyReal8ToPtr(o2,o2_pr,nx1*nx2*nx3*nx4*nx5*nx6*nx7)
    call mxCopyReal8ToPtr(o3,o3_pr,nx1*nx2*nx3*nx4*nx5*nx6*nx7)
    
    ! Deallocate arrays
    deallocate(x8)
    deallocate(x9)
    deallocate(pf1)
    deallocate(pf2)
    deallocate(pf3)
    deallocate(o1) 
    deallocate(o2) 
    deallocate(o3)    
end subroutine mexFunction
subroutine allterp(	nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9, &
					x8,x9, &
					x8i,x9i, &
					pf1,pf2,pf3, &
					o1,o2,o3)
    implicit none
    mwSize :: nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9
    double precision :: x8i,x9i,x8(nx8),x9(nx9)
    double precision, dimension(nx1,nx2,nx3,nx4,nx5,nx6,nx7,nx8,nx9) :: pf1,pf2,pf3
    double precision, dimension(nx1,nx2,nx3,nx4,nx5,nx6,nx7) :: o1,o2,o3
    
    double precision :: s8, s9
    double precision :: x8i_min, x9i_min
    mwSize loc8, loc9
    double precision, dimension(2) :: xi, xi_left, xi_right, w_2, w_1
	double precision, dimension(2) :: w8, w9
	mwSize :: m8, m9
	double precision :: wtemp
	s8 = x8(2) - x8(1)
	s9 = x9(2) - x9(1)	
	x8i_min = x8i - x8(1)
	loc8 = min(nx8-1,max(1,floor(x8i_min/s8) + 1))
	
	x9i_min = x9i - x9(1)
	loc9 = min(nx9-1,max(1,floor(x9i_min/s9) + 1))
	xi = [x8i, x9i]
	xi_left = [x8(loc8), x9(loc9)]
	xi_right = [x8(loc8+1), x9(loc9+1)]
	w_2 = (xi - xi_left)/(xi_right - xi_left)
	w_1 = 1 - w_2	
	w8 = [w_1(1), w_2(1)]
	w9 = [w_1(2), w_2(2)]
	
	o1 = 0.d0
	o2 = 0.d0
	o3 = 0.d0	
	do m9 = 0,1
    	    do m8 = 0,1
			  wtemp = w8(m8+1)*w9(m9+1)
			  o1 = o1 + wtemp*pf1(:,:,:,:,:,:,:,loc8+m8,loc9+m9)
			  o2 = o2 + wtemp*pf2(:,:,:,:,:,:,:,loc8+m8,loc9+m9)
			  o3 = o3 + wtemp*pf3(:,:,:,:,:,:,:,loc8+m8,loc9+m9)
            end do
        end do			  
end subroutine allterp
The code is a part of an optimization routine to interpolate policy functions. This is the first time I have tried to compile a mex file. Any help would be greatly appreciated.

Best Answer

Three comments:
1) You should NEVER, NEVER, NEVER pass literal integers into mex API functions. It can often be the case that your compiler will treat those literal integers as INTEGER*4 type, when in fact the MATLAB version & API library you are linking to is expecting INTEGER*8 type. Arguments are not type promoted in Fortran like they are in C/C++. To avoid this potential problem, always used typed variables to pass arguments into mex API functions using the exact signatures specified in the doc. So this
plhs(1) = mxCreateNumericArray(4,[nx1,nx2,nx3,nx4,nx5,nx6,nx7],myclassid,0)
should be something like this instead
    mwSize ndim
    integer*4 ComplexFlag
        :
    ndim = 4
    ComplexFlag = 0
    plhs(1) = mxCreateNumericArray(ndim,[nx1,nx2,nx3,nx4,nx5,nx6,nx7],myclassid,ComplexFlag)
That being said, the above code doesn't make any sense to me as written. Why are you passing in a 7-element dimension array and telling mxCreateNumericArray that you are only passing in 4 dimensions? Those last three dimensions of nx5, nx6, nx7 will be ignored as the code is currently written, possibly creating a situation where the result is undersized and downstream code will try to access invalid memory and crash. Seems like that first argument should be a 7, not a 4.
2) There are no checks in the code to see if the dimension values passed in actually match the data variables passed in. If there is a mismatch you can easily access invalid memory and crash MATLAB. That is, there are no checks like this in the code:
mwPointer, external :: mxGetNumberOfElements
    :
if( mxGetNumberOfElements(prhs(14)) /= nodes ) then
    call mexErrMsgTxt("Number of elements mismatch for prhs(14)")
endif
if( mxGetNumberOfElements(prhs(15)) /= nodes ) then
    call mexErrMsgTxt("Number of elements mismatch for prhs(15)")
endif
if( mxGetNumberOfElements(prhs(16)) /= nodes ) then
    call mexErrMsgTxt("Number of elements mismatch for prhs(16)")
endif
3) You don't have basic checks for input integrity in your code for number and type of inputs. E.g., you should have stuff like this up front (not necessarily an inclusive list):
integer k
    :
if( nlhs /= 3 ) then
    call mexErrMsgTxt("Need to specify three outputs")
endif
if( nrhs /= 16 ) then
    call mexErrMsgTxt("Need 16 full double real inputs")
endif
do k=1,16
    if( .not.mxIsDouble(prhs(k)) .or. mxIsSparse(prhs(k)) .or. mxIsComplex(prhs(k)) ) then
        call mexErrMsgTxt("All inputs must be full double real")
    endif
enddo
Etc.
If the types and sizes and sparsity and complexity etc. are not what the mex routine expects, this can easily crash MATLAB. You should make all of these checks before you access any variables and before you allocate any variables.
Related Question