MATLAB: Fortran MEX returns inconsistent answers with -largeArrayDims

fortranmex

I programmed a MEX function that gives reasonable results most of the time, but about once every 5 runs or so it will return a matrix of NaNs, or a matrix of huge numbers (1.0e+43 when it should be around 100 or so). When I rerun the function with the exact same inputs the problem often goes away. I'm guessing that it has to do with MEX interfacing. Does anyone have any idea what could be going on? Here is the MEX interface:

#define __LP64__
#include "fintrf.h"
        SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
        IMPLICIT NONE
        MWPOINTER :: PLHS(*),PRHS(*)
        INTEGER*4 :: NLHS,NRHS ! REMAINS THE SAME FOR 64BIT system
  !-----------------------------------------------------------------------
        MWSIZE,PARAMETER::dp=SELECTED_REAL_KIND(12, 60)
        MWSIZE :: ndim, err
        MWSIZE :: dims(6)
        MWSIZE :: num_choice_params, num_HC_params, size_Emax_poly
        MWSIZE :: num_ages, num_J, num_cohorts, num_K, num_T, num_Q, num_HC_shocks
        MWSIZE, EXTERNAL :: MXGETM, MXGETN
        MWPOINTER ::Choice_Pr, HC_Pr, Wage_Pr, EMAX_Pr
        MWPOINTER, EXTERNAL :: MXGETPR,MXCREATENUMERICARRAY
        REAL(dp), ALLOCATABLE :: Choice_Params(:), HC_params(:), Wage_stacked(:,:), EMAX(:,:,:,:,:, :)
  ! Arguments for mxCreateNumericArray
        INTEGER*4 :: CLASSID, COMPLEXFLAG, MXCLASSIDFROMCLASSNAME 
! ASSIGN POINTERS TO THE VARIOUS PARAMETERS 
      Choice_pr = MXGETPR(PRHS(1)) 
      HC_pr = MXGETPR(PRHS(2)) 
      Wage_pr = MXGETPR(PRHS(3)) 
      size_Emax_poly = MXGETM(PRHS(4)) 
      num_J = MXGETM(PRHS(3)) 
      num_K = MXGETM(PRHS(5)) 
      num_Q = MXGETM(PRHS(6))   
      num_HC_shocks= MXGETM(PRHS(7))   
      num_choice_params = MXGETM(PRHS(1)) 
      num_HC_params = MXGETM(PRHS(2))   
     !   n = MXGETN(PRHS(1))
       num_ages=66-18+1 
       num_T=2012-1994+1
       num_cohorts=num_ages+num_T-1
! ALLOCATE SIZES TO INPUT AND OUTPUT MATRICES
        ALLOCATE( Choice_Params(num_choice_params), STAT = err )
        IF (err .NE. 0) THEN
          CALL MEXERRMSGTXT('MultMexError: Out of memory Choice params')
        ENDIF
        ALLOCATE( HC_Params(num_HC_params), STAT = err )
        IF (err .NE. 0) THEN
          CALL MEXERRMSGTXT('MultMexError: Out of memory HC_params')
        ENDIF
        ALLOCATE( Wage_stacked(num_J, num_T*num_K), STAT = err )
        IF (err .NE. 0) THEN
          CALL MEXERRMSGTXT('MultMexError: Out of memory Wage_stacked')
        ENDIF
        ALLOCATE( EMAX(num_ages, num_J, num_cohorts, num_Q, num_K+1, 3*size_Emax_poly+1), STAT = err)
        IF (ERR .NE. 0) THEN
            CALL MEXERRMSGTXT('MultMexError: Out of memory EMAX')
        END IF
! COPY RIGHT HAND ARGUMENTS TO LOCAL ARRAYS
        CALL MXCOPYPTRTOREAL8(Choice_pr, Choice_params, num_choice_params)
        CALL MXCOPYPTRTOREAL8(HC_pr, HC_params, num_HC_params) 
        CALL MXCOPYPTRTOREAL8(Wage_pr, Wage_stacked, num_t*num_j*num_k)         
! OUTPUT MATRIX AND POINTER
        CLASSID = MXCLASSIDFROMCLASSNAME('double')
        COMPLEXFLAG = 0
        ndim = 6
        dims(1) = num_ages
        dims(2) = num_J
        dims(3) = num_cohorts
        dims(4) = num_Q
        dims(5) = num_K+1
        dims(6) = 3*size_Emax_poly+1
        PLHS(1) = MXCREATENUMERICARRAY(ndim,dims,CLASSID,COMPLEXFLAG)
        Emax_pr = MXGETPR(PLHS(1))
! CREATE 3D MATRIX 
        CALL Calc_Emax(Choice_params,HC_Params,Wage_stacked, size_Emax_poly, num_HC_shocks, &
         num_ages, num_J, num_cohorts, num_K, EMAX, num_choice_params, num_HC_params, num_T, num_Q )
! COPY RESULTS TO OUTPUT POINTER 
        CALL MXCOPYREAL8TOPTR(EMAX,EMAX_pr, num_ages*num_J* num_cohorts* (num_K+1)*(3*size_Emax_poly+1)*num_Q)
        DEALLOCATE(Choice_params)
        DEALLOCATE(HC_params)
        DEALLOCATE(Wage_stacked)
        RETURN
        END SUBROUTINE MEXFUNCTION

Best Answer

I haven't looked at the code in detail, but one thing I did notice was the declarations for the MXGETM and MXGETN functions are not up to date. It should be using MWPOINTER, not MWSIZE. E.g.,

MWPOINTER, EXTERNAL :: MXGETM, MXGETN

See the doc:

http://www.mathworks.com/help/matlab/apiref/mxgetm.html

http://www.mathworks.com/help/matlab/apiref/mxgetn.html

Related Solutions

MATLAB: How does one return allocated memory from a Fortran Mex function

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.

MATLAB: Building Mex file with Fortan

Your main problem is that you need to create the output array PLHS(1) before you access it. Also your NLHS and NRHS types are wrong. Try this:

#include <fintrf.h>
C
C Gateway routine
C
    SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
    implicit none
C mexFunction arguments:
    MWPOINTER PLHS(*), PRHS(*)
    INTEGER NLHS, NRHS
    mwPointer, external :: mxGetPr, mxGetPi
    mwPointer, external :: mxCreateDoubleMatrix, mxCreateNumericArray
    mwSize, external :: mxGetM, mxGetN
    MWPOINTER YPP 
    mwSize M, N
    integer*4 ComplexFlag
    REAL *8 RYPP(1)
    RYPP(1) = 10
    M = 1
    N = 1
    ComplexFlag = 0
    PLHS(1) = mxCreateDoubleMatrix( M, N, ComplexFlag )
    YPP = MXGETPR(PLHS(1))
    CALL MXCOPYREAL8TOPTR(RYPP, YPP, M*N)
    RETURN
    END

Best Answer

Related Solutions

MATLAB: How does one return allocated memory from a Fortran Mex function

MATLAB: Building Mex file with Fortan

Related Question