Inverse matrix in fortran mexFunction - matlab lapack library

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

[Muthu Annamalai]

Have you run it on a debugger and gotten a stack trace? It seems to me like you are having a memory issue here.

Answer 1

[James Tursa]

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.

[Tue]

Thaks for answer. This solved the issue.

Answer 2

[zohar]

Hi don't know anything with fortran.

1) From my experience, every time I used LAPACK I checked info.

2) What is the type of your matrix (general, symmetric positive-definite..), I see it complex.

3) And the big question WHY DO YOU WANT TO FIND INV ?? I am sure you do not need to do that ! Use Factorize matrix like ZGETRF and then Solve system like ZGETRS.

4) ZGETRF and ZGETRS depend on the type of your matrix.

5) Check this line: plhs(1) = mxCreateDoubleMatrix(m,n,1). m != n !?!?

[Tue]

Hi zohar. To clarify: 1) I should. 2) The matrix to invert is general complex (nonsymmetric). 3) I need the inverse because the resulting matrix is needed in > 20 different matrix products afterwards. 4) I will look this up. 5) n=m always.

[Muthu Annamalai]

Zohar means that your matrix maybe ill-conditioned, i.e. your det(A) \approxly 0, and inverse could be troublesome. Using LU, factoriazation and other associated transform methods are more reliable.