Crash during mxCopyReal8toPtr in Fortan 90 Mex File when Arrays are Large

1 view (last 30 days)
My code crashes when the inputs ( math_input1 and math_input2 ) are large arrays (in other words, when ntimes is large), resulting in large output arrays. It crashes during the calls to mxCopyReal8toPtr. When I comment out those four calls, the code completes without error (though obviously the mxArrays have only been instantiated and the math results have not been copied over). I clearly have enough memory to perform the computation and create the mxArrays large enough to contain the results, but I am guessing that the mxCopy* might require significantly more memory in order to copy the results over to the Matlab arrays? If this is the case and it's just a memory issue, I want to understand the requirements of mxCopy* in order to prevent code crashes. It would be even better if there was an alternative that didn't have the same limitations.
I have stripped most of the functionality and the actual computation out of my code, but I have maintained the same structure, variable declarations, variable allocations, and matlab API calls. I have also confirmed that this stripped code exhibits the same behavior as my actual routines.
Run and compiled on a 64-bit machine, called using Matlab R2011b.
GATEWAY ROUTINE BELOW
#include "fintrf.h"
subroutine MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
use module_stripped
implicit none
mwPointer PLHS(*), PRHS(*), mxGetData, mxGetM, mxGetN, pr
mwPointer mxCreateDoubleMatrix, mxCreateNumericArray
mwSize mxGetNumberOfElements, NLHS, NRHS, i, mxClassIDFromClassName, nfreqs, nevents, ntimes, neventbins, nfreqsxtimes, neventsxeventbins
mwSize, parameter :: one = 1
mwSize, allocatable :: m(:),n(:),e(:)
integer(4) :: int4zero
parameter (int4zero = 0)
real(8) :: ti, dt
real(8), allocatable, dimension(:) :: math_input1, math_input2, math_input1_full, math_input2_full
real(8), allocatable, dimension(:) :: math_result1, math_result2, math_result3, math_result4
real(8), allocatable, dimension(:) :: freqs
real(8), allocatable, dimension(:,:) :: t1, t2
allocate(m(NRHS),n(NRHS),e(NRHS))
do i = one, NRHS
m(i) = mxGetM(PRHS(i))
n(i) = mxGetN(PRHS(i))
e(i) = mxGetNumberOfElements(PRHS(i))
end do
nfreqs = e(10)
nevents = m(7)
neventbins = n(7)
ntimes = e(4)
nfreqsxtimes = nfreqs*ntimes
neventsxeventbins = nevents*neventbins
allocate(math_input1_full(ntimes))
allocate(math_input2_full(ntimes))
allocate(t1(nevents,neventbins))
allocate(t2(nevents,neventbins))
call mxCopyPtrToReal8(mxGetData(PRHS(4)),math_input1_full,ntimes)
call mxCopyPtrToReal8(mxGetData(PRHS(14)),math_input2_full,ntimes)
call mxCopyPtrToReal8(mxGetData(PRHS(5)),dt,one)
call mxCopyPtrToReal8(mxGetData(PRHS(6)),ti,one)
call mxCopyPtrToReal8(mxGetData(PRHS(7)),t1,neventsxeventbins)
call mxCopyPtrToReal8(mxGetData(PRHS(8)),t2,neventsxeventbins)
if ( nint((maxval(t2(1:nevents,1:neventbins))-ti)/dt+1.0)<ntimes ) then
ntimes = nint( (maxval(t2(1:nevents,1:neventbins))- ti)/dt+1.0)
end if
allocate(math_input1(ntimes),math_input2(ntimes))
math_input1(1:ntimes) = math_input1_full(1:ntimes)
math_input2(1:ntimes) = math_input2_full(1:ntimes)
deallocate(math_input2_full,math_input1_full)
allocate(math_result1(nfreqsxtimes))
allocate(math_result2(nfreqsxtimes))
allocate(math_result3(nfreqsxtimes))
allocate(math_result4(nfreqsxtimes))
allocate(freqs(nfreqs))
! 2
! 3
! 4
! 5
! 6
! 7
! 8
! 9
call mxCopyPtrToReal8(mxGetData(PRHS(10)),freqs,nfreqs)
! 11
! 12
! 13
PLHS(1) = mxCreateDoubleMatrix(ntimes,nfreqs,int4zero)
PLHS(2) = mxCreateDoubleMatrix(ntimes,nfreqs,int4zero)
PLHS(3) = mxCreateDoubleMatrix(ntimes,nfreqs,int4zero)
PLHS(4) = mxCreateDoubleMatrix(ntimes,nfreqs,int4zero)
call module_subroutine(math_result1, math_result2, math_result3, math_result4, math_input2, math_input1, dt, ti, t1, t2, nfreqs, freqs, nevents, ntimes)
pr = mxGetData(PLHS(1))
call mxCopyReal8ToPtr(math_result1,pr,nfreqsxtimes)
deallocate(math_result1)
call mxCopyReal8ToPtr(math_result2,mxGetData(PLHS(2)),nfreqsxtimes)
deallocate(math_result2)
call mxCopyReal8ToPtr(math_result3,mxGetData(PLHS(3)),nfreqsxtimes)
deallocate(math_result3)
call mxCopyReal8ToPtr(math_result4,mxGetData(PLHS(4)),nfreqsxtimes)
deallocate(math_result4)
return
end subroutine
COMPUTATIONAL ROUTINE BELOW IN A SEPARATE FILE
#include "fintrf.h"
module module_stripped
implicit none ! no types based on naming
public module_subroutine
private module_lower_subroutine
contains
RECURSIVE subroutine module_subroutine(math_result1, math_result2, math_result3, math_result4, math_input2, math_input1, dt, ti, t1, t2, nfreqs, freqs, nevents, ntimes)
implicit none
integer ERR
character*128 MSG
!integer(4) :: nfreqs, nf, nj, ntimes, nevents
mwSize :: nfreqs, nf, nj, ntimes, nevents
real(8), intent(in) :: dt, ti, math_input1(*), math_input2(*)
real(8) :: math_result1(*), math_result2(*), math_result3(*), math_result4(*)
real(8), allocatable, dimension(:) :: temp1, temp2, temp3, temp4
real(8) freqs(*), t1(*), t2(*)
!$OMP PARALLEL PRIVATE(temp1,temp2,temp3,temp4)
if (.NOT. allocated(temp1))then
allocate(temp1(ntimes),STAT=ERR,ERRMSG=MSG)
end if
if (.NOT. allocated(temp2))then
allocate(temp2(ntimes),STAT=ERR,ERRMSG=MSG)
end if
if (.NOT. allocated(temp3))then
allocate(temp3(ntimes),STAT=ERR,ERRMSG=MSG)
end if
if (.NOT. allocated(temp4))then
allocate(temp4(ntimes),STAT=ERR,ERRMSG=MSG)
end if
!$OMP DO
do nf = 1,nfreqs
call module_lower_subroutine(temp1,temp2,temp3,temp4, math_input2, math_input1, dt, ti, ntimes)
math_result1((nf-1)*ntimes+1:(nf-1)*ntimes+ntimes) = temp1(1:ntimes)
math_result2((nf-1)*ntimes+1:(nf-1)*ntimes+ntimes) = temp2(1:ntimes)
math_result3((nf-1)*ntimes+1:(nf-1)*ntimes+ntimes) = temp3(1:ntimes)
math_result4((nf-1)*ntimes+1:(nf-1)*ntimes+ntimes) = temp4(1:ntimes)
end do
!$OMP END DO
if (allocated(temp1))then
deallocate(temp1)
end if
if (allocated(temp2))then
deallocate(temp2)
end if
if (allocated(temp3))then
deallocate(temp3)
end if
if (allocated(temp4))then
deallocate(temp4)
end if
!$OMP END PARALLEL
end subroutine
RECURSIVE subroutine module_lower_subroutine(math_result1, math_result2, math_result3, math_result4, math_input2, math_input1, dt, ti, ntimes)
use omp_lib
implicit none
!integer(4) :: ntimes
mwSize :: ntimes
real(8) :: math_input2(*), math_input1(*), ti, dt
real(8) :: math_result1(*), math_result2(*), math_result3(*), math_result4(*)
math_result4(1:ntimes) = 0.0
math_result3(1:ntimes) = 0.0
math_result2(1:ntimes) = 0.0
math_result1(1:ntimes) = 0.0
end subroutine
end module module_stripped

Answers (1)

James Tursa
James Tursa on 4 Mar 2014
Edited: James Tursa on 4 Mar 2014
I haven't looked at your code in any great detail, but one thing I noticed was this:
call mxCopyPtrToReal8(mxGetData(PRHS(5)),dt,1)
call mxCopyPtrToReal8(mxGetData(PRHS(6)),ti,1)
You should never pass constant literals to a MATLAB API function in Fortran, because the type can change between 32-bit and 64-bit versions. You can get away with this in C/C++ when the argument is passed by value (a copy is made and automatically converted to the correct type if necessary per the prototype), but you can't get away with this in Fortran where arguments like the above 1 are typically passed by reference. In the above lines, the signature from the doc is this:
subroutine mxCopyPtrToReal8(px, y, n)
mwPointer px
real*8 y(n)
mwSize n
In particular, the last argument is mwSize, which is typically a 32-bit integer on 32-bit systems and a 64-bit integer on 64-bit systems. But constant literal integers such as 1 are typically 32-bit integers on both 32-bit systems and 64-bit systems. If that happens to be true on your system, you would be passing addresses to 32-bit integers to a routine that is expecting addresses to 64-bit integers. Of course, subsequent code will then access invalid memory and then lots of bad things can happen (garbage integer value being used, crash, etc.). Try this instead:
mwSize, parameter :: one = 1
:
call mxCopyPtrToReal8(mxGetData(PRHS(5)),dt,one)
call mxCopyPtrToReal8(mxGetData(PRHS(6)),ti,one)
  4 Comments
James Tursa
James Tursa on 5 Mar 2014
I will take another look. One other signature mismatch to be aware of:
integer nlhs, nrhs
You have them as mwSize.
James Tursa
James Tursa on 5 Mar 2014
Edited: James Tursa on 5 Mar 2014
It might be a good time to start putting in more argument checking up front as well. E.g.,
if( mxIsDouble(PRHS(4)) /= 1 .OR. mxIsDouble(PRHS(5)) /= 1 .OR. etc )
call mexErrMsgTxt("Inputs are not double class as required")
endif
Also put in checks to ensure the number of elements you are copying are in fact present in the input arrays you are copying from and that there is enough space in the target as well.

Sign in to comment.

Categories

Find more on Fortran with MATLAB in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!