MATLAB: MEX file (FORTRAN) calling a LAPACK routine crashed

mex fortran lapack

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

Best Answer

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.
Related Question