Segmentation violation error when I run a mex file

I need to call a Fortran file when I run matlab. The fortran file was provided by roms.org. I can use mex to generate .mexa64 using

mex -f ./mexopts.sh -v mexrect.F

However, when I run it in matlab by typing

z = mexrect (1i, 65, 2, 33, 34,65);

a very serious error appears, showing "I need to call a Fortran file when I run matlab. The fortran file was provided by roms.org. I can use mex to generate mexrect.mexa64 using

mex -f ./mexopts.sh -v mexrect.F

However, when I run it in matlab by typing

z = mexrect (1i, 65, 2, 33, 34,65);

a very serious error appears, showing "Segmentation violation detected at Wed Aug 7 21:21:26 2013" and then matlab crashes.

Does anyone know why?

Thanks.

Li

Below is the Fortran file (mexrect.F).

#include <fintrf.h>
c
c mex routine implementing rect.f as taken from Ives, D.C. and
c R. M. Zacharias "Conformal mapping and Orthogonal Grid
c Generation", aiaa-87-2057.
c
c
c 
c The calling sequence from MATLAB should be
c
c   >> zn = mexrect ( z, n, n1, n2, n3, n4 );
c
c So 
c    1)  nlhs = 1
c    2)  plhs = 
c    3)  nrhs = 6
c    4)  prhs(1) = z
c        prhs(2) = n
c        prhs(3) = n1
c        prhs(4) = n2
c        prhs(5) = n3
c        prhs(6) = n4
c

      subroutine mexFunction ( nlhs, plhs, nrhs, prhs )

      Implicit Double Precision (A-H,O-Z)
  c        implicit none

      MWPOINTER plhs(*), prhs(*)
  c      integer*4 plhs(*), prhs(*)
      integer*4 nlhs, nrhs
  c     integer*8 nlhs, nrhs

c
c     Size of input arguments
    integer size_m, size_n

c
c     Grid outline to be transformed.
    Complex*16 z(40000)

c
c     MWPOINTERs to real and imaginary parts of z input argument,
c     which is prhs(1)
    MWPOINTER zpr, zpi
c      integer*4 zpr,zpi
c
c     MWPOINTERs to real parts of n, n[1234] input arguments, which
c     are prhs(2) thru prhs(6)
    MWPOINTER np, n1p, n2p, n3p, n4p
c      integer*4 np, n1p, n2p, n3p, n4p
c
c     Size of input z in terms of matlab dimensions, i.e.,
c     #rows x #columns
c      integer zsize
    mwSize zsize
    mwSize numel
    mwPointer, external :: mxGetM, mxGetN
      mwSize zsize
    mwSize numel
    mwPointer, external :: mxGetM, mxGetN, mxGetPi, mxGetPr
    mwPointer, external :: mxIsComplex, mxIsNumeric
    mwPointer, external :: mxClassIdFromClassName
    mwPointer, external :: mxCreateNumericMatrix700
c 
c     These point to the real and imaginary parts of the 
c     matlab structure plhs(1)
    MWPOINTER zlpr, zlpi
c      integer*4 zlpr,zlpi

c     We apparently have to copy the matlab arguments to
c     real arrays first, then cast them to integers.
    real*8 nr(1), n1r(1), n2r(1), n3r(1), n4r(1)
    integer n, n1, n2, n3, n4

c     Check for proper number of arguments.
    if ( nrhs .ne. 6 ) then
       call mexErrMsgTxt('mexrect requires 6 arguments')
    elseif ( nlhs .ne. 1 ) then
       call mexErrMsgTxt ( 'Only one output argument allowed' )
    endif

c     Check to see that the first input was complex.
c     This is the "Z" argument.
    if ( mxIsComplex(prhs(1)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( 'First argument to mexrect must be complex' )
    endif

c     Check to see that the 2nd thru 6th arguments are numeric.
c     This is "N", "N1", "N2", "N3", and "N4".
    if ( mxIsNumeric(prhs(2)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( '2nd argument to mexrect must be numeric' )
    elseif ( mxIsNumeric(prhs(3)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( '3rd argument to mexrect must be numeric' )
    elseif ( mxIsNumeric(prhs(4)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( '4th argument to mexrect must be numeric' )
    elseif ( mxIsNumeric(prhs(5)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( '5th argument to mexrect must be numeric' )
    elseif ( mxIsNumeric(prhs(6)) .ne. 1 ) then
       call mexErrMsgTxt 
   +        ( '6th argument to mexrect must be numeric' )
    endif

c
c     Check to see that the 2nd thru 6th arguments are scalars.
    size_m = mxGetM ( prhs(2) )
    size_n = mxGetN ( prhs(2) )
    if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
       call mexErrMsgTxt
   +        ('2nd argument to mexrect must be a scalar' )
    endif

      size_m = mxGetM ( prhs(3) )
      size_n = mxGetN ( prhs(3) )
    if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('3rd argument to mexrect must be a scalar' )
      endif

      size_m = mxGetM ( prhs(4) )
      size_n = mxGetN ( prhs(4) )
    if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('4th argument to mexrect must be a scalar' )
      endif

      size_m = mxGetM ( prhs(5) )
      size_n = mxGetN ( prhs(5) )
    if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('5th argument to mexrect must be a scalar' )
      endif

      size_m = mxGetM ( prhs(6) )
      size_n = mxGetN ( prhs(6) )
    if ( size_n .ne. 1 .or. size_m .ne. 1 ) then
         call mexErrMsgTxt
     +        ('6th argument to mexrect must be a scalar' )
      endif

c
c     Create matrix for the return argument.
    size_m = mxGetM ( prhs(1) )
    size_n = mxGetN ( prhs(1) )
    zsize = size_m*size_n
    plhs(1) = mxCreateNumericMatrix(size_m, size_n,
   +                                mxClassIdFromClassName('double'),
   +                                1 );
    zlpr = mxGetPr(plhs(1))
    zlpi = mxGetPi(plhs(1))

c
c     Load the data into Fortran matrices.
    zpr = mxGetPr(prhs(1))
    zpi = mxGetPi(prhs(1))
    call mxCopyPtrToComplex16(zpr,zpi,z,zsize)

      numel=1;

      np = mxGetPr(prhs(2))
    call mxCopyPtrToReal8(np,nr,numel)
    n = nr(1)

      n1p = mxGetPr(prhs(3))
    call mxCopyPtrToReal8(n1p,n1r,numel)
    n1 = n1r(1)

      n2p = mxGetPr(prhs(4))
    call mxCopyPtrToReal8(n2p,n2r,numel)
    n2 = n2r(1)

      n3p = mxGetPr(prhs(5))
    call mxCopyPtrToReal8(n3p,n3r,numel)
    n3 = n3r(1)

      n4p = mxGetPr(prhs(6))
    call mxCopyPtrToReal8(n4p,n4r,numel)
    n4 = n4r(1)

c
c     Call the computational subroutine.
    call rect ( z, n, n1, n2, n3, n4 )

c
c     Load the output into a MATLAB array.
    call mxCopyComplex16ToPtr(z,zlpr,zlpi,zsize)
    return
    end

c
    Subroutine RECT(Z,N,N1,N2,N3,N4)
c
c     this subroutine is taken from ives,d.c. and
c     r.m.zacharias "conformal mapping and orthogonal grid generation"
c     aiaa-87-2057.
c
c     The only modification is the addition of the "tracking_error"
c     warning message.  This allows matlab to gracefully exit a call
c     to RECT that is going badly.
c
    Implicit Double Precision (A-H,O-Z)
    Complex*16 Z(1), Z0, ZD
    Double Precision R(40000), T(40000)
    PI = 4.D0 * DATAN(1.D0)
    Do 70 I = 1, N
      IM = N - MOD(N-I+1,N)
      IP = 1 + MOD(I,N)
      ZD = (Z(IM)-Z(I)) / (Z(IP)-Z(I))
      ALPHA = DATAN2(DIMAG(ZD),DREAL(ZD))
      If (ALPHA.LT.0) ALPHA = ALPHA + PI + PI
      PWR = PI / ALPHA
      If (I.EQ.N1.OR.I.EQ.N2.OR.I.EQ.N3.OR.I.EQ.N4) Then
        ZD = (Z(IM)-Z(I)) / CDABS(Z(IM)-Z(I))
        Do 10 J = 1, N
          Z(J) = DCMPLX(0.D0,1.D0) * Z(J) / ZD
 10     Continue
        PWR = PWR / 2.
      End If
      H1IN = 100.
      H1AX = -100.
      TP = 0.
      Do 40 J = 2, N
        ZD = Z(MOD(J+I-2,N)+1) - Z(I)
        R(J) = CDABS(ZD)
        T(J) = DATAN2(DIMAG(ZD),DREAL(ZD)) - PI - PI - PI - PI - PI -
   *        PI
        Do 20 K = 1, 7
          If (DABS(T(J)-TP).LT.PI) Go To 30
          T(J) = T(J) + PI + PI
 20     Continue
       call mexErrMsgTxt ( 'warning - tracking error ' )
     return
c        pause ' warning - tracking error '
c        Call CRASH(TTYOUT,4,0)
 30     TP = T(J)
        H1AX = DMAX1(H1AX,T(J)*PWR)
        H1IN = DMIN1(H1IN,T(J)*PWR)
 40   Continue
      PWR = DMIN1(PWR,1.98D0*PI*PWR/(H1AX-H1IN))
      Z(I) = DCMPLX(0.D0,0.D0)
      Do 50 J = 2, N
        Z(MOD(J+I-2,N)+1) = R(J) ** PWR * 
   *        CDEXP(DCMPLX(0.D0,T(J)*PWR))
 50   Continue
      ZD = 1. / (Z(N2)-Z(N1))
      Z0 = Z(N1)
      Do 60 J = 1, N
        Z(J) = (Z(J)-Z0) * ZD
 60   Continue
 70 Continue
    Return
    End

[Jan]

The initial part of the message seems to be duplicated accidentally.

Without a proper formatting the code is not readable. So please follow the "? Help" link to learn, how to present your question in a form, which allows to understand it.

Answer 1

[James Tursa]

You should use "implicit none" in your code ... it is one of the few friends you have in Fortran.

You should declare all of your functions to make sure the compiler knows the type that is returned.

You should never use constants for arguments to the API functions ... always use variables to ensure that the argument types match.

E.g.,

The documentation shows the following signature for mexFunction:

subroutine mexFunction(nlhs, plhs, nrhs, prhs)
integer*4 nlhs, nrhs
mwPointer plhs(*), prhs(*)

So use integer*4 for the nlhs and nrhs, not integer. Probably not a source of the error since I would expect integer to be integer*4 even on a 64-bit system, but why not use the published signature?

Then consider mxGetM and mxGetN. Since you don't declare them they will get an integer return type by default (probably integer*4 on your system), but it is entirely possible that the actual return type is integer*8 on your 64-bit system. So declare them properly to avoid this potential system crashing mistake:

mwPointer, external :: mxGetM, mxGetN

Then consider the mxCopyPtrToReal8 calls you make using a constant in the last argument. That constant will likely be an integer*4, but the signature is:

subroutine mxCopyPtrToReal8(px, y, n)
mwPointer px
real*8 y(n)
mwSize n

That mwSize might be an integer*8 on a 64-bit system. So this call:

call mxCopyPtrToReal8(np,nr,1)

should look like this instead:

mwSize numel
    :
numel = 1
call mxCopyPtrToReal8(np,nr,numel)

I can't comment on anything else you have posted until you format it correctly with the "{ } Code" button. It is unreadable until you do this.

[Jan]

[EDITED, Relocated from answer section, Jan] Qiang wrote:

Hi James,

Many thanks for your advices.

After I declared all the mx* function, it works.

I added below in the code:

      mwSize zsize
      mwSize numel
      mwPointer, external :: mxGetM, mxGetN, mxGetPi, mxGetPr
      mwPointer, external :: mxIsComplex, mxIsNumeric
      mwPointer, external :: mxClassIdFromClassName
      mwPointer, external :: mxCreateNumericMatrix700

Cheers,

Li

[Jan]

@Qiang: I've moved your comment to the comment section. Now you can select another accepted answer.

Answer 2

[Qiang]

Hi James,

I am dealing with mex a lot recently. The previous code works with your advices. However, another code doesn't work. Do you have some time to take a look?

Thanks,

Li

Below is the .F code. There are also some extra subroutines. I think they do not matter, so I didn't copy here. If you need them, I will copy them to you.

#include <fintrf.h>

           subroutine mexfunction(nlhs, plhs, nrhs, prhs)
       implicit none
C-----------------------------------------------------------------------
C     (integer) Replace integer by integer*8 on the DEC Alpha and the
C     SGI 64-bit platforms
C
      mwPointer plhs(*), prhs(*)
      integer*8 x,y,z,datain,xhat,yhat,zhat
      integer*8 pmap,dataout,errout,a,b
      mwPointer, external :: mxGetPr, mxCreateDoubleMatrix
C-----------------------------------------------------------------------
C
! matlab call:
!     [dataout,errout]=oaTest(x,y,data,xhat,yhat);
      integer*8 nlhs, nrhs
      mwPointer, external :: mxGetM, mxGetN
      mwSize m_in, n_in, size
      mwSize PtIn,PtOut,PtIns, PtZin,PtZout       

! get address of input      arrays

        x = mxGetPr(prhs(1))
        y = mxGetPr(prhs(2))
    z = mxGetPr(prhs(3))
    datain = mxGetPr(prhs(4))               
        xhat = mxGetPr(prhs(5))
        yhat = mxGetPr(prhs(6))            
        zhat = mxGetPr(prhs(7))                
    pmap = mxGetPr(prhs(8))    
        a = mxGetPr(prhs(9))
        b = mxGetPr(prhs(10))
  
  !set output matrix according to sizes of Sn (3)

      m_in = mxGetM(prhs(4))
      n_in = mxGetN(prhs(4))    
      PtIn = m_in 
  PtZin= n_in
      m_in = mxGetM(prhs(7))
      n_in = mxGetN(prhs(7))    
  PtOut = m_in 
  PtZout= n_in
      plhs(1) = mxCreateDoubleMatrix(m_in, n_in, 0)
      dataout = mxGetPr(plhs(1))  
  n_in = 1
      plhs(2) = mxCreateDoubleMatrix(m_in, n_in, 0)  
  errout = mxGetPr(plhs(2))

      m_in = mxGetM(prhs(8))
      n_in = mxGetN(prhs(8))    
  PtIns = n_in

      call oa_ex(%val(x), %val(y),%val(z),%val(datain),
     c               %val(xhat),%val(yhat),%val(zhat),
     c           PtIn,PtOut,PtZin,PtZout,%val(dataout),
     c               %val(errout),%val(a), %val(b),
     c               PtIns, %val(pmap))

       return
   end