Problem passing large array between matlab and mex function

[EDIT: 20110523 19:35 CDT - reformat - WDR]

On 64-bit Windows 7 and 64-bit matlab, I have the following problem. When I pass a 2D array of type single (float32) and size (43500, 29200) from matlab to a mex function, only the first 196458176 values are passed to the mex function correctly. The same problem occurred when I passed a 2D array of type single and size (43500, 29200) from the mex function to matlab. The mex function is written in Fortran and the code was compiled by the Intel 64 Fortran compiler XE version 12.0.4.196 Build 20110427.

I can send you the Fortran code (io_64bit.f) for the mex function. I test the mex function by the following script:

nC = 43500;
nR = 29200;
din = ones(nC, nR, 'single');
[dout1, dout2] = io_64bit(din);
sdin   = sum(sum( int32(din) ));
sdout1 = sum(sum( int32(dout1) ));
sdout2 = sum(sum( int32(dout2) ));

If all goes well, sdout1 = sdin = - sdout2 = 1270200000 (= nC * nR)

Following is the Fortran code:

c***********************************************************************
c
c     Check fidelity of large array passed from Matlab to mex function
c     and vice versa.
c
c     [dout1, dout2] = io_64bit(din)
c
c***********************************************************************
c
c     Original: 18 May 2011                Final: 18 May 2011
c
c***********************************************************************
      subroutine mexFunction(nlhs, plhs, nrhs, prhs)
      implicit none

        integer*4    nlhs, nrhs
  c     <-----------------------------------------------------------------
  c     For 32-bit OS, change following to integer*4
  c     --------------------------------------------
        integer*8    plhs(*), prhs(*)
        integer*8    mxGetCell, mxGetM, mxGetN, mxGetPr
        integer*8    mxGetNumberOfElements, mxClassIDFromClassName
        integer*8    mxCreateNumericArray, mxCreateNumericMatrix
        integer*8    mxIsInt8, mxIsInt32, mxIsSingle
        integer*8    mexPrintf

        integer*8    nCol, nRow, din_p, dout1_p, dout2_p

        integer*8    n, ClassID, ComplexFlag
  c     ----------------------------------------------------------------->

        character (len = 48) eMsg
        integer    i, iStatus
        real       xmin, xmax

        real,      allocatable :: din(:,:)

        logical, parameter :: lDebug = .true.

c     ------------------------------------------------------------------
      if(nrhs /= 1) call mexerrmsgtxt('1 input argument is required')
      if(nlhs /= 2) call mexerrmsgtxt('2 outputs are produced')

        nCol = mxGetM(prhs(1))
        nRow = mxGetN(prhs(1))

        n = mxGetNumberOfElements(prhs(1))

        if(lDebug) then
           open(1,file='mex_io.log',status='unknown')
           write(1,*) nCol, nRow, n
        endif

        if(n /= nCol*nRow) then
           call mexErrMsgTxt('Cannot determine size of 2D array dem')
        endif

c     -------------------------------- Check inputs are of correct types
      eMsg = 'Argument x must be of type single / float32'
      do i=1, 1
         if(mxIsSingle(prhs(i)) == 0) then
            write(eMsg(10:10),'(i1)') i
            call mexErrMsgTxt(eMsg)
         endif
      enddo

c     -------------------------------- Allocate memory and setup pointer
      allocate(din(nCol,nRow),stat=iStatus)
      call check_allocation(iStatus, 'din')

        din = 0.0

        if(lDebug) write(1,*) 'Done memory allocation'

        din_p    = mxGetPr(prhs(1))

        if(lDebug) write(1,*) 'Done pointer ', din_p

c     -------------------------------------------- Copy data from Matlab
      call mxCopyPtrToReal4(din_p, din, n)

        if(lDebug) then
           write(1,*) 'Done copy ptr to real'
           call minmax(n, din, xmin, xmax)
           write(1,*) xmin, xmax
        endif

c     ------------------------------ Create output matrix of type single
      ClassID = mxClassIDFromClassName('single')
      ComplexFlag = 0

        plhs(1) = mxCreateNumericMatrix(nCol, nRow, ClassID, ComplexFlag)
        plhs(2) = mxCreateNumericMatrix(nCol, nRow, ClassID, ComplexFlag)

        dout1_p  = mxGetPr(plhs(1))
        dout2_p  = mxGetPr(plhs(2))

        if(lDebug) then
           write(1,*) ClassID
           close(1)
        endif

c     -------------------------------------------- Copy data 1 to Matlab
      call mxCopyReal4ToPtr(din, dout1_p, n)

c     -------------------------------------------- Copy data 2 to Matlab
      din = -1.0
      call mxCopyReal4ToPtr(din, dout2_p, n)

        deallocate(din)

        return
        end subroutine mexFunction
  c
  c=======================================================================
  c
        subroutine minmax(n, x, xmin, xmax)
        integer*8 n
        real xmin, xmax, x(n)

        integer i, ix, m

        xmin = huge(xmin)
        xmax = -xmin

        m = 0
        do i=1, n
           xmin = min(xmin, x(i))
           xmax = max(xmax, x(i))
           if(x(i) /= 0.0) ix = i
           if(x(i) == 0.0) m = m + 1
        enddo

c     ------------------------------------------------------------------
c     ix = no of values passed from Matlab correctly
c
c     m  = no of no-values from Matlab
c
c     (n - m) should = ix
c
c     ------------------------------------------------------------------
      write(1,*) ix, m, n - m

        return
        end subroutine minmax
  c
  c=======================================================================
  c
        subroutine check_allocation(iStatus, name)
        integer        iStatus
        character*(*)  name

        if(iStatus /= 0) then
           call mexErrMsgTxt('Cannot allocate memory for array '//name)
           stop
        endif

        return
        end subroutine check_allocation

[Walter Roberson]

You said "if all goes well", but what output do you get instead?

Answer 1

[Philip Borghesani]

Your FORTRAN code is missing:

#include "fintrf.h"

At the beginning of the file. This header is needed when using -largeArrayDims because it it responsible for switching which mx* functions are called. See the standard mex documentation and this solution be sure to read the faq that is linked to at the end of solution. The last section of the text faq has useful FORTRAN information.

[James Tursa]

Yes, I know ... sorry about that. It was pointed out some time ago and appears in the comments sections of download page. I have been working on updated files for some time (fixes that problem and also problems related to pointing to 0-sized arrays and the LOC function not being implemented properly on some compilers) but have yet to upload it. In particular, there will likely be some problems with the examples on 64-bit systems. I hope to get the new code uploaded soon.

[Philip Borghesani]

James I was not commenting on your code, I believe the missing include is Joseph's problem.

[James Tursa]

Ah! Thanks for the clarification. (Well, if you ever do download my code you now know about some problems that need fixing)

[Joseph]

Thanks Philip and James. Now I get a better understanding of the mex gateway between Matlab and your own computer codes. Also, I was initially puzzled by "-largeArrayDims flag does not change the way your code is compiled". The statement is true as is. What it does not say represents the crux of problem I am having. The inclusion of fintrf.h and pre-processing your Fortran code by fpp before it is compiled changes your code to call '64-bit' subroutines. I tend to compile mex function outside Matlab and I also tend not to rely fpp. This is purely a personal choice.

Answer 2

[Titus Edelhofer]

Hi,

did you compile the code using the "-largeArrayDims" flag? See also http://www.mathworks.com/support/solutions/en/data/1-5C27B9/?solution=1-5C27B9

Titus

[Joseph]

I did try compiling the code using the -largeArrayDims flag and without using the flag. In both cases, the results are the same.

sdout1 = - sdout2 = 196458176

Please note that my array size of 43500 * 29200 < 2^32-1.

I was also told by Matlab support that using the -largeArrayDims flag does not change the way the code is compiled!

[Philip Borghesani]

-largeArrayDims does not change how the code is compiled but it does change which mx* functions are called and how they must be called. To work with large arrays the flag is needed.

[Joseph]

Hi Philip, that is something I do not understand how this works. If you mex code is not compiled and linked with the correct mx* functions, how does the mexw64 that you compiled know which function to call?

Answer 3

[Walter Roberson]

Your array has fewer than 2^32-1 elements, but because each element is 4 bytes, your array has more than 2^32 bytes. The excess above 2^32 corresponds exactly to 196458176 array elements.

[Joseph]

Thanks Walter for pointing this out. May be the documentation on this should point this out as you described.

Answer 4

[James Tursa]

The first thing I would do would be to examine all of the function/subroutine interfaces that you use and make sure they are correct. e.g., in the online doc there is this:

integer*4 mxClassIDFromClassName(classname)
character*(*) classname

So this function returns an integer*4 type, yet in your code you have it declared as returning an integer*8 type. Same comment for mxIsSingle, mxIsInt8, and mxIsInt32. So double check the doc for your particular installation to make sure you have all of the function interfaces declared properly.

The next thing I would ask is why you are not using the mwPointer and mwSize types that MATLAB gives you for this. e.g., I would typically do something like this:

#ifndef mwSize
#define mwSize integer*4
#endif

#ifndef mwPointer
#define mwPointer integer*4
#endif
     :
     :
mwPointer plhs(*), prhs(*)
mwPointer, external :: mxGetCell
mwPointer, external :: mxGetPr
mwSize, external :: mxGetM
mwSize, external :: mxGetN
mwSize, external :: mxGetNumberOfElements
integer*4, external :: mxClassIDFromClassName
mwPointer, external :: mxCreateNumericArray
mwPointer, external :: mxCreateNumericMatrix
integer*4, external :: mxIsInt8
integer*4, external :: mxIsInt32
integer*4, external :: mxIsSingle
integer*4, external :: mexPrintf

mwSize nCol, nRow
mwPointer din_p, dout1_p, dout2_p
mwSize n
integer*4 ClassID, ComplexFlag

The last thing I would point out is that you can avoid all of your huge data copying with the use of the %val construct in your function calls. This will be a much better use of resources and will run quite a bit faster as well. e.g.,

!      allocate(din(nCol,nRow),stat=iStatus)
!      call check_allocation(iStatus, 'din')
!        din = 0.0
!       if(lDebug) write(1,*) 'Done memory allocation'
      din_p    = mxGetPr(prhs(1))
      if(lDebug) write(1,*) 'Done pointer ', din_p
!     -------------------------------------------- Copy data from Matlab
!      call mxCopyPtrToReal4(din_p, din, n)
        if(lDebug) then
!           write(1,*) 'Done copy ptr to real'
           call minmax(n, %val(din_p), xmin, xmax)
           write(1,*) xmin, xmax
        endif

How much RAM do you have?

[Joseph]

Hi James, Thanks for your suggestions. On 64-bit Matlab, mx functions need to be declared integer*8 (see, eg, yprimefg.F).

I will use %val construct. Thanks for that.

I have 48 GB of RAM.

[James Tursa]

I would have expected all of the mx functions that return pointers to be integer*8 on a 64-bit system. I was more wondering about the others that should have been integer*4 & that perhaps they could be messing you up. Are you saying that your specific doc says that mxClassIDFromClassName (and others) returns an integer*8 on your system? That would be news to me, and would necessitate me changing my Fortran 95 MATLAB interface code on the FEX to accommodate it. Normally I would have suggested you use my FEX code for the Fortran interface, but unfortunately it only supports double (real*8) variables at present.

[James Tursa]

P.S. There *is* a way with the Intel compiler to get Fortran pointers to mxArray data areas without using %VAL like in my FEX submission. It basically relies on the fact that a Fortran pointer to a scalar is represented by the Intel compiler as just a simple address (no special descriptor) and makes an implicit call to a routine that relies on that representation to return a 2D Fortran pointer to the same data area. Let me know if you are interested and I can post the code.