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.