Fortran MEX returns inconsistent answers with -largeArrayDims

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

Answer 1

[James Tursa]

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

[mycolas]

James,

Thanks for the response. I made the change and I'm still getting the same issue. Is it possible there's a problem with how I'm declaring my inputs as double precision? Currently for the argument choice_params I have:

MWSIZE,PARAMETER::dp=SELECTED_REAL_KIND(12, 60)

MWPOINTER ::Choice_Pr

REAL(dp), ALLOCATABLE :: Choice_Params(:)

Choice_pr = MXGETPR(PRHS(1))

CALL MXCOPYPTRTOREAL8(Choice_pr, Choice_params, num_choice_params)

Choice_params is then declared as dp within the fotran subroutine.

It seems like the parameters are sometimes off by a factor of 1000- so I expect a value of 1.76 but instead I get back 1760. Maybe there's confusion caused by my "dp" and MXCOPYPTRTOREAL8 expecting something else?

[James Tursa]

Since you said your results were correct most of the time, I didn't think dp was the problem. That is, if you are passing in double data and dp is actually single inside your mex routine, I wouldn't have expected your results to be correct at all. But I don't know what Calc_Emax is doing so maybe that is not a valid observation on my part. So maybe you could just force the issue and hard code a REAL(8) in your code instead of using your REAL(dp) and see what happens.

[James Tursa]

Another thing that might be tripping you up is this expression as one of your arguments:

num_ages*num_J* num_cohorts* (num_K+1)*(3*size_Emax_poly+1)*num_Q

If the default integer size for literal integers is larger than MWSIZE, then this expression will result in a larger integer than MWSIZE being passed in. I always prefer to force the type in these cases so I don't have to worry about that. I.e., create another variable to use specifically for the argument. E.g.,

MWSIZE :: num
    :
num = num_ages*num_J* num_cohorts* (num_K+1)*(3*size_Emax_poly+1)*num_Q)
CALL MXCOPYREAL8TOPTR(EMAX,EMAX_pr, num)

Also, double check the array sizes in your Calc_Emax routine to make sure they match what is being passed in. Can you post the Calc_Emax signature stuff? Or maybe the entire routine if it isn't too long?

[mycolas]

James, I figured it out. I hadn't fully assigned one of my arrays in Calc_Emax and so it was spitting out wonky answers.

Thank you so much for your help, I really appreciate it. Working through it with you lead me right to the mistake. If I graduate on time you'll definitely be a reason why!