Fortran MEX file memory problem

2013 2 28
1 个回答
回答已采纳
9 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

James Tursa
James Tursa 2013-2-28

0 个投票

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

3 个评论
显示 1更早的评论 隐藏 1更早的评论

Giorgio
Giorgio 2013-2-28
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
James Tursa 2013-2-28
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
Giorgio 2013-3-1
编辑：Giorgio 2013-3-1
I follow your hint, and now it works fine.
Thank you for your help!
Giorgio

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Write C Functions Callable from MATLAB (MEX Files) 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by