Fortran MEX file memory problem

Dear all,

I implemented an iterative solver, to solve the generic complex linear system A*x = b, by using a Fortran subroutine, while my main program is all written in Matlab (actually I'm using Matlab R2010b). Now, when my A matrix exceeds a certain dimension (about 1500x900 double complex values) I receive this message from the Matlab editor:

Caught MathWorks::System::FatalException

I made also a fortran implementation of my main and, after increasing the stack size (under properties --> Linker --> System --> Stack Reserve size) on my Microsoft Visual Studio 2008 everything works fine (with any dimension of A matrix). But I can't do this when I generate my mex file by using Matlab.

In the following the Gateway routine which call my fortran subroutine:

if true #include "fintrf.h"
!======================================================================
#if 0
!     
!     z_gatewayGMRES.f
!     .f file needs to be preprocessed to generate .for equivalent
!     
#endif
!     
!======================================================================
!
!     Gateway routine
    subroutine mexFunction(nlhs, plhs, nrhs, prhs)
!     Declarations
    implicit none
!
!----------------------------------------------------------------------
!     mexFunction arguments:
!     (pointer) Replace integer by integer*8 on 64-bit platforms
    integer*8   PLHS(*)
    integer*8   PRHS(*)
!
    integer*8   nlhs
    integer*8   nrhs
!
!----------------------------------------------------------------------
!     Function declarations:
!     (pointer) Replace integer by integer*8 on 64-bit platforms
    integer*8   mxCreateDoubleMatrix
    integer*8   mxCreateNumericMatrix
    integer*8   mxGetData
    integer*8   mxGetPr, mxGetPi
    integer*8 mxGetM, mxGetN
    integer*8 mxClassIDFromClassName
!
!----------------------------------------------------------------------
!  I/O & LOCAL VARIABLE DEFINITIONS
! ---------------------------------------------------------------------
!
    INTEGER*8   mA, nA, sizeA, mb
    INTEGER*8   maxoutPin, tolPin, omgPin
    INTEGER*8   ApIn, bPin, nbinPin
! Input --------------------------------------------------------------- C
    COMPLEX*16, ALLOCATABLE ::   A(:,:), b(:)
    INTEGER*8   nbin, maxout
    REAL*8      nbinR, maxoutR
    REAL*8      tol, omg
    REAL*8      t_in, t_out
! Output -------------------------------------------------------------- C
    INTEGER*8   conv, nbout
    COMPLEX*16, ALLOCATABLE ::   x(:)
    REAL*8,     ALLOCATABLE ::   res(:)
! Local Variables
    INTEGER*8   classid
    INTEGER*8   complexflag
    INTEGER*8   r_datasize, z_datasize
    INTEGER*8   m, n
!
! --------------------------------------------------------------------- C
!
#if defined MSWIND
! For Windows only!
! This resets the floating point exception to allow divide by zero,
! overflow and invalid numbers. 
!
    INTEGER(2) CONTROL
    CALL GETCONTROLFPQQ(CONTROL)
    CONTROL = CONTROL .OR. FPCW$ZERODIVIDE
    CONTROL = CONTROL .OR. FPCW$INVALID
    CONTROL = CONTROL .OR. FPCW$OVERFLOW
    CALL SETCONTROLFPQQ(CONTROL)
#endif
!
! CHECK THE DIMENSIONS
!
    mA    = mxGetM(PRHS(1))
    nA    = mxGetN(PRHS(1))
    sizeA = mA*nA
!
    mb    = mxGetM(PRHS(2))
!
    IF( mA .NE. mb ) THEN
       CALL mexErrMsgTxt('The number of rows is not the same')
    ELSE
       m = mA
       n = nA
    ENDIF
!
! ---------------------------------------------------------------------------- 
!
    r_datasize  = 8
    z_datasize  = 16
!
! ASSIGN POINTERS TO THE INPUT PARAMETERS
!
    nbinPin     = mxGetData(PRHS(3))
    maxoutPin   = mxGetData(PRHS(4))
    tolPin      = mxGetData(PRHS(5))
    omgPin      = mxGetData(PRHS(6))
!
! INITIALIZE HEAP MEMORY FOR INPUT MATRIX/ARRAYS
!
    ALLOCATE(A(m,n))
    ALLOCATE(b(m))
!
! ----------------------------------------------------------------------------
! COPY RIGHT HAND ARGUMENTS TO LOCAL ARRAYS OR VARIABLES
!
    CALL mxCopyPtrToReal8(nbinPin, nbinR, 1)
    CALL mxCopyPtrToReal8(maxoutPin, maxoutR, 1)
    nbin     = INT(nbinR)
    maxout   = INT(maxoutR)
    CALL mxCopyPtrToReal8(tolPin, tol, 1)
    CALL mxCopyPtrToReal8(omgPin, omg, 1)
!
    CALL mxCopyPtrToComplex16( mxGetPr(prhs(1)), &
                               mxGetPi(prhs(1)),A,sizeA)
    CALL mxCopyPtrToComplex16( mxGetPr(prhs(2)), &
                               mxGetPi(prhs(2)),b,m)
!
! ----------------------------------------------------------------------------
!     INITIALIZE HEAP MEMORY FOR LEFT HAND ARGUMENTS
!
    ALLOCATE(x(n))
    ALLOCATE(res(maxout))
! ----------------------------------------------------------------------------
!  Iterative solver
!
    CALL It_solver(m, n,                &
                   A, b, tol,           &
                   nbin, maxout, omg,   &
                   x, conv, res, nbout, &
                   t_in, t_out          )
!
! ----------------------------------------------------------------------------
!     CREATE RETURN ARGUMENT ARRAYS
!
    classid  = mxClassIDFromClassName('int64')
    complexflag = 1
    PLHS(1) = mxCreateDoubleMatrix(n, 1, complexflag)
    complexflag = 0
    PLHS(2) = mxCreateNumericMatrix(1, 1, classid, complexflag)
    PLHS(3) = mxCreateDoubleMatrix(maxout, 1, complexflag)
    PLHS(4) = mxCreateNumericMatrix(1, 1, classid, complexflag)
!
!     Load the output into a MATLAB array.
!
    CALL mxCopyComplex16ToPtr(x,  mxGetPr(PLHS(1)), &
                                  mxGetPi(PLHS(1)),n)
    CALL mxCopyInteger4ToPtr(conv, mxGetData(PLHS(2)), 1)
    CALL mxCopyReal8ToPtr(res,mxGetPr(PLHS(3)),maxout)
!
! ---------------------------------------------------------------------------- 
!
    DEALLOCATE(A)
    DEALLOCATE(b)
    DEALLOCATE(x)
    DEALLOCATE(res)
!
! ----------------------------------------------------------------------------
!
    RETURN
    END
  % code
end

Any suggestion will be really appreciated.

Thank you very much in advance,

Giorgio

Answer 1

[James Tursa]

Can you post the variable declarations for the subroutine It_solver? Is there a large variable there that gets allocated off of the stack? I.e., a local variable that gets its dimension from m or n?

Also, all of those mxCopyEtc calls make my brain hurt. You might take a look at this Fortran 95 API interface submission on the FEX:

http://www.mathworks.com/matlabcentral/fileexchange/25934-fortran-95-interface-to-matlab-api-with-extras

[Giorgio]

Thank you for your fast answer; in the following the variable declaration for the subroutine It_solver:

if true
    subroutine It_solver(m, n,                &
                         A, b, tol,           &
                         nbin, maxout, omg,   &
                         x, conv, res, nbout, &
                         t_in, t_out          )
!
! ====================================================================
! Arguments
!
!  Input:
!        m:             number of rows (i.e. the number of equations)
!        n:             number of columns (i.e. the nubmer of unknowns)
!        A(m,n):        input matrix
!        b(1:m):        right hand side
!        tol:           tolerance to stop the algorithm
!        nbin:          number of inner iterations
!        maxout:        maximum number of outter iterations
!        omg:           acceleration parameter for the SOR method
!
!  Output:
!        x(1:n):        the solution of the linear least squares problem
!        conv:          the status of the computation, see below
!        res(1:maxout): norm of the residual r=A^(T)b-A^(T)Ax
!        nbout:         total number of outter iterations necessary
!                       to find the solution vector
! =========================================================================
!
    USE ConstantsModule
!
    IMPLICIT NONE
!
! Input -----------------------------------------------------------------
!
    integer(kind = 4),intent(in)  :: nbin, maxout, m, n
    real(kind = 8),   intent(in)  :: tol, omg      
    complex(kind = 8),intent(in)  :: A(m,n), b(m)
!
! Output ----------------------------------------------------------------
!
    integer(kind = 4),intent(out) :: conv, nbout
    real(kind = 8),   intent(out) :: res(maxout)
    real(kind = 8),   intent(out) :: t_in, t_out
    complex(kind = 8),intent(out) :: x(n)
!
! Work arrays -----------------------------------------------------------
!
    integer(kind = 4) :: i, j, k, l
    integer(kind = 4) :: nbin_opt
    real(kind = 8)    :: omg_opt
    real(kind = 8)    :: beta, r_inprod, r_tmp, loc_tol
    real(kind = 8)    :: t_in1, t_in2, t_out1, t_out2
    real(kind = 8)    :: Aei(n), c(maxout)
    complex(kind = 8) :: d, z_inprod, z_tmp
    complex(kind = 8) :: w(n),y(maxout),r(m),s(maxout),g(maxout+1)
    complex(kind = 8) :: H(maxout+1,maxout), V(n,maxout+1)
!========================================================================
!
    ...
!
!  ----------------------------------------------------------------------
!
    end subroutine It_solver
end

I also try to use the mxCalloc() function and the %val() to pass variable to the subroutine, but the result is still the same.

Thank you for your kind help.

Giorgio

[James Tursa]

In the subroutine It_solver, the work arrays Aei(n), c(maxout), etc all get their memory off of the stack. I would suggest making all of them allocatable instead since that will result in them getting their memory from the heap (much larger than the stack). E.g., turn this:

real(kind = 8) :: H(maxout+1,maxout)

into this

real(kind = 8), allocatable :: H(:,:)
    :
allocate(H(maxout+1,maxout))
    :
deallocate(H)

... and so on for all the other local dimensioned work arrays.

[Giorgio]

I follow your hint, and now it works fine.

Thank you for your help!

Giorgio