MEX file (FORTRAN) calling a LAPACK routine crashed

I am now trying to learn some Fortran MEX coding using MATLAB 2010B on 64-bit Linux. I start with the simplest case, i.e. MATLAB calling a LAPACK subroutine, zlarfg. The code test.f90 is shown as follws. But I find it is very strange that the compilation command

       mex -largeArrayDims -lifcore test.f90 -lmwlapack -lmwblas

results in a mex file which never works and makes Matlab GUI crash.

      SUBROUTINE MEXFUNCTION( nlhs, plhs, nrhs, prhs )

      IMPLICIT NONE
!
!     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
      COMPLEX*16        CZERO
      PARAMETER         ( CZERO = ( 0.0D0, 0.0D0 ) )
!
!     .. Mex-file interface parameters ..
      integer nlhs, nrhs
      mwPointer plhs(*), prhs(*)
      mwSignedIndex   mxGetM, mxGetN, M, N
      mwPointer   mxCreateDoubleMatrix, mxGetPr, mxGetPi
      integer   mxIsNumeric, mxIsComplex
      integer           LDA, IP, istate
!
!     .. Allocatable arrays ..
      complex*16, Allocatable :: AZ(:,:), TAU(:)
      COMPLEX*16        ALPHA, tau1

      IP = 1
      N = mxGetN( PRHS(IP) )
      M = mxGetM( PRHS(IP) )
      IF ( mxIsNumeric( PRHS(IP) ) == 0 )  CALL mexErrMsgTxt( 'A must be a numerical matrix' )
      IF ( M /= N )    CALL mexErrMsgTxt( 'A must be a square matrix' )
      IF ( mxIsComplex( PRHS(IP) ) == 0 )  CALL mexErrMsgTxt( 'A must be a complex matrix' )

      LDA =  N

 ALLOCATE ( AZ(LDA,LDA), TAU(LDA), STAT = istate)

      CALL mxCopyPtrToComplex16( mxGetPr(PRHS(1)), mxGetPi(PRHS(1)), AZ, N*N)
      TAU = CZERO

      IP = N
      ALPHA = AZ( IP - 1, IP )

      call zlarfg(IP-1, ALPHA, AZ(1:(IP-2), IP), 1, tau1)

TAU(IP-1) = tau1
open(unit=90, file='tmp1.txt')

write(90,*) 'AZ',AZ(1:IP, IP)
WRITE(90,*) 'TAU', TAU1
close(90)

      PLHS(1) = mxCreateDoubleMatrix( N, N, 1 )
      CALL mxCopyComplex16ToPtr(AZ, mxGetPr(PLHS(1)),mxGetPi(PLHS(1)), N*N)

     PLHS(2) = mxCreateDoubleMatrix( N, 1, 1 )
     CALL mxCopyComplex16ToPtr(TAU, mxGetPr(PLHS(2)),mxGetPi(PLHS(2)), LDA)

DEALLOCATE (AZ, TAU)

      END

Answer 1

[James Tursa]

Fortran does not automatically convert arguments to the correct type like C/C++ does. You need to get the arguments the same as the signature requires. In the case of 64-bit BLAS and LAPACK libraries for MATLAB, that typically means using mwSignedIndex for all integer arguments. So in this call:

    call zlarfg(IP-1, ALPHA, AZ(1:(IP-2), IP), 1, tau1)

The IP-1 and the 1 are likely 32-bit integers, even on a 64-bit machine, while zlarfg is likely expecting 64-bit integer arguments. Change these to mwSignedIndex.

Also, you've got LDA and 1's in some of the MATLAB API calls. Again, you may be passing 32-bit integers to routines expecting 64-bit integers. It is critical in Fortran to use the exact types in the signatures for API calls ... in your case you should be using mwSize variables. mxGetM and mxGetN also are mistyped (you have mwSignedIndex while the signature has mwPointer), but you will probably get away with that one since they are likely the same.

Bottom line: Scrub your code to make sure all routine calls are using the correct argument types.

Answer 2

[Tian]

Thanks very much for your answer. Following your suggestion, the resulting mex file works properly.

I am very sorry to say that due to lack of knowledge about the syntax of mexfunction, I try to copy other's example on the internet, so I can not figure out what is the wrong with my code, even if sometimes the code happened to run properly.

BTW, what does mwSignedIndex, mwSignedIndex and mwPointer exactly mean? And where should I use them when writing a Fortran mexfunction?

[James Tursa]

The simple advice is always look at the doc for the signature of the function you are calling, and use the exact variable types in the call. For Fortran in particular, the mwPointer type is an integer that is the same size as a pointer, and is used to hold a pointer value. So mwPointer will be a 32-bit integer in a 32-bit app and it will be a 64-bit integer in a 64-bit app. Thus it is convenient to use in the API to mimic a size_t in C for an integer in certain function calls (even though the name suggests pointer).

[Tian]

Thanks again for your suggestion. However, just now, I tried a slightly more complicated case, finding that your answer alone seemingly could not make the following code run properly. Specifically, the code breaks down when calling ZDOTU, at which I again felt quite confused.

 #include "fintrf.h"
      SUBROUTINE MEXFUNCTION( nlhs, plhs, nrhs, prhs )

      IMPLICIT NONE

      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
      COMPLEX*16        CZERO
      PARAMETER         ( CZERO = ( 0.0D0, 0.0D0 ) )
      mwSize, parameter :: ione = 1

      integer nlhs, nrhs
      mwPointer  plhs(*), prhs(*)
      mwPointer  mxGetM, mxGetN, mxCreateDoubleMatrix, mxGetPr, mxGetPi
      integer   mxIsNumeric, mxIsComplex

      mwSize      M, N, LDA, LDWORK, LZWORK, ILO, IP
      integer     istate
      CHARACTER         JOBQ

      REAL*8, Allocatable :: DWORK(:)
      complex*16, Allocatable :: AZ(:,:), D(:), ZWORK(:), W(:,:)

      IF (NRHS /= ione) CALL mexErrMsgTxt('requires only one input' )

      IP = ione
      N = mxGetN( PRHS(IP) )
      M = mxGetM( PRHS(IP) )
      IF ( mxIsNumeric( PRHS(IP) ) == 0 )  CALL mexErrMsgTxt( 'A must be a numerical matrix' )
      IF ( M /= N )    CALL mexErrMsgTxt( 'A must be a square matrix' )
      IF ( mxIsComplex( PRHS(IP) ) == 0 )  CALL mexErrMsgTxt( 'A must be a complex matrix' )

      LDA =  N
      !

         LDWORK = N - ione
         LZWORK =  N - ione

open(unit=100, file='tmp.txt')
write(100, *) M
write(100, *) N
write(100, *) LDWORK
write(100, *) LZWORK
close(100)

ALLOCATE (W(LDA,LDA), ZWORK(LZWORK), AZ(LDA,LDA), DWORK(LDWORK), stat = istate)

 open(unit=120, file='tmpaz.txt', position = 'append')
 write(120,*) 'istate', istate

      ZWORK = CZERO
      DWORK = 0.d0
      W = CZERO
      AZ = CZERO

       CALL mxCopyPtrToComplex16( mxGetPr(PRHS(ione)), mxGetPi(PRHS(ione)), AZ, N*N)

write(120,*) 'AZ', AZ
close(120)

       call syax('U', LDA, LDA, AZ, LDA, DWORK, ZWORK, W, LDA)

open(unit=90, file='tmp1.txt')
write(90,*) 'AZ', AZ
close(90)

      PLHS(ione) = mxCreateDoubleMatrix( N, N, ione)
      CALL mxCopyComplex16ToPtr(AZ, mxGetPr( PLHS(ione) ), &
                                      mxGetPi( PLHS(ione) ), N*N )

   DEALLOCATE( W, AZ, ZWORK, DWORK)

      END

      SUBROUTINE SYAX( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
      IMPLICIT NONE
      mwSize , parameter::  ione =  1
      COMPLEX*16,parameter:: CZERO=( 0.0D0, 0.0D0 ), &
                             ONE= ( 1.0D0, 0.0D0 )

      CHARACTER          UPLO
      mwSize::            LDA, LDW, N, NB

      COMPLEX*16         A( LDA, N ), TAU( N - ione ), W( LDW, NB)
      DOUBLE PRECISION   E( N-ione )

      mwSize ::           I, IW
      COMPLEX*16         ALPHA

      EXTERNAL           ZAXPY, ZGEMV, ZSYMV, ZLACGV, ZLARFG, ZSCAL

      LOGICAL            LSAME
      COMPLEX*16         ZDOTU
      EXTERNAL           LSAME, ZDOTU

      open(unit = 100,file = 'inax.txt')
       write(100,*) 'uplo',uplo
        write(100,*) 'NB', NB

    E = 0.d0,  
       W = czero,
       I = N

       ALPHA = A( I-ione, I)

        CALL ZLARFG(I-ione, ALPHA, A(ione, I), ione, TAU(I-ione))

         write(100,*) 'ALPHA', ALPHA
         write(100,*) 'TAU(I-1)',TAU(I-ione)

         IW = N

         CALL ZSYMV( 'Upper', I-ione, ONE, A, LDA, A( ione, I ), ione, &
                           ZERO, W( ione, IW ), ione )

         CALL ZSCAL( I-ione, TAU( I-ione ), W( ione, IW ), ione )

          write(100,*) 'W(1,IW)',W(ione,IW)
          ALPHA = ZDOTU(I-ione, A(ione, I), ione, W(ione, IW), ione)
          write(100,*) 'ALPHA', ALPHA
           close(100)

              END

[James Tursa]

You are still passing mwSize variables (I, IONE) instead of mwSignedIndex variables. Try fixing that first.