MATLAB: Inverse matrix in fortran mexFunction – matlab lapack library

fortranlapackmex

I am trying to use LAPACK to invert a matrix in Fortran through mex. Depending on the matrix dimension Matlab crashes. E.g. for Hin=diag(1:3);[T1]=TestInvMex(int32(1),Hin) gives correct answer. Hin=diag(1:6);[T1]=TestInvMex(int32(1),Hin) crashes. Matlab version 2012b on linux. My code TestInvMex.F90:
#include "fintrf.h"
!======================================================================
! mex -v -largeArrayDims TestInvMex.F90 -output TestInvMex FOPTIMFLAGS='-O3' -lmwlapack -lmwblas
!======================================================================
!     Gateway routine
      subroutine mexFunction(nlhs, plhs, nrhs, prhs)
!     Declarations
      implicit none
!     mexFunction arguments:
      mwPointer plhs(*), prhs(*)
      integer nlhs, nrhs
!     Function declarations:
      mwPointer mxCreateDoubleMatrix, mxGetPr, mxGetPi
      integer mxIsNumeric
!      mwPointer mxGetM, mxGetN, mxIsInt8, mxIsInt32
      mwSignedIndex mxGetM, mxGetN, mxIsInt8, mxIsInt32
!     Pointers to input/output mxArrays:
      mwPointer mode_pr, H_RE_pr, H_IM_pr, Hinv_RE_pr, Hinv_IM_pr
!     Array information:
!      mwsize m, n
      mwSignedIndex m, n
!     Define variables - Arguments for computational routine:
      integer, parameter :: dp=SELECTED_REAL_KIND(15) ! dp=15 digits, sp->6 instead.
      integer, parameter :: I32=SELECTED_INT_KIND(15) ! INT64=15, INT32=6.
      mwSignedIndex,parameter :: myOne=1,myTwo=2,myFour=4
      INTEGER(I32) :: mode 
      COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: H
      COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Hinv
      COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: work  ! work array for LAPACK
      integer :: info
      integer, DIMENSION(:), ALLOCATABLE :: ipiv   ! pivot indices
!     Get the size of the input arrays (m rows, n collums).
      m = mxGetM(prhs(2))
      n = mxGetN(prhs(2))
!     Allocate:
      ALLOCATE(H(m,n),Hinv(m,n),work(m),ipiv(m))
!     Create matrix for the return argument. NB:1 to have complex elements
      plhs(1) = mxCreateDoubleMatrix(m,n,1)
      Hinv_RE_pr = mxGetPr(plhs(1))
      Hinv_IM_pr = mxGetPi(plhs(1))
!     Input: 
      mode_pr = mxGetPr(prhs(1))
      H_RE_pr = mxGetPr(prhs(2))
      H_IM_pr = mxGetPi(prhs(2))
!     Load the data into Fortran arrays.
      call mxCopyPtrToComplex16(H_RE_pr,H_IM_pr,H,m*n)
      call mxCopyPtrToInteger4(mode_pr,mode,myOne)
!     Compute inverse:
      Hinv=H ! prevent H from being overwritten by LAPACK
      call ZGETRF(n, n, Hinv, n, ipiv, info) ! Complex dp
      call ZGETRI(n, Hinv, n, ipiv, work, n, info) ! Complex dp
!     Load the data into the output to MATLAB.
      call mxCopyComplex16ToPtr(Hinv,Hinv_RE_pr,Hinv_IM_pr,m*n)
        return
        end

Best Answer

The general crashing problem is caused by the fact that you use complex LAPACK routines when you don't have complex inputs or outputs. So you are copying stuff into and out of NULL pointers. E.g.,
      plhs(1) = mxCreateDoubleMatrix(m,n,1)
      Hinv_RE_pr = mxGetPr(plhs(1))
      Hinv_IM_pr = mxGetPi(plhs(1))
Thus Hinv_IM_pr will be a NULL pointer (i.e., no memory allocated for imaginary part).
Then downstream you copy into the Hinv_IM_pr pointer, which is invalid:
      call mxCopyComplex16ToPtr(Hinv,Hinv_RE_pr,Hinv_IM_pr,m*n)
Similar situation applies to your input. E.g.,
      H_IM_pr = mxGetPi(prhs(2))
If the input is real, as in your example, H_IM_pr will be a NULL pointer (invalid, no memory allocated for imaginary part). Then you try to copy from it here:
      call mxCopyPtrToComplex16(H_RE_pr,H_IM_pr,H,m*n)
In both cases cited above, I would expect MATLAB to crash because you are accessing an invalid pointer.
Also I would question how you are setting up your integer arguments. E.g., you have an int32(1) as an input argument, then you do this in your code to recover it:
      integer, parameter :: I32=SELECTED_INT_KIND(15) ! INT64=15, INT32=6.
      mwSignedIndex,parameter :: myOne=1,myTwo=2,myFour=4
      INTEGER(I32) :: mode 
          :
      mwSignedIndex,parameter :: myOne=1,myTwo=2,myFour=4
          :
      mode_pr = mxGetPr(prhs(1))
          :
      call mxCopyPtrToInteger4(mode_pr,mode,myOne)
It seems to me that you are using machine specific code here in an attempt to be robust, but it looks rather convoluted to me. The whole process of getting a value into an integer scalar from MATLAB can be done in one line:
      mode = mxGetScalar(prhs(1))
The above will work for any numeric class of input, not just an int32.
Also, info and ipiv should be mwSignedIndex, not integer, since they are being passed into the LAPACK routine.
As a final general comment, it is advised to put in a lot of class and dimension and complex checks into your code up front to make sure the input is as expected before copying from the data pointers. That way you could handle both real and complex inputs in the same routine by selectively copying from the imaginary part only if it is actually physically present.
Related Question