segfault in matrix inversion using LAPACK in mex

Question

I'm trying to do a matrix inversion in a FORTRAN code which is called in MATLAB. But when I run the code in debug mode it can call and compute "dgetrf" without any problem but then I get a segfault in "dgetri". Can anyone give me an insight about any possible reason for such an error?

Here is the gateway function for the test function I use:

#include "fintrf.h"

!
SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
use mexf90
implicit none
integer*8, intent(in)   :: PRHS(*)      ! Pointers carrying the input data
integer*8, intent(out)  :: PLHS(*)      ! Pointers carrying the output data
integer*4, intent(in)   :: NLHS,NRHS    ! REMAINS THE SAME FOR 64BIT system
!-----------------------------------------------------------------------
integer*8               :: err, nA, mA
integer*8, pointer      :: matAR, AInvR
character(200)          :: errMsg
integer*4               :: txt, classId
integer*4, external     :: mexprintf

! THE PART BELOW SEPARATED WITH "!$$" IS FOR WINDOWS ONLY
! THIS RESETS THE FLOATING POINT EXCEPTION
! TO ALLOW DIVIDE BY ZERO,
! OVERFLOW AND INVALID 

!$$
!#if defined MSWIND
!        INTEGER(2) CONTROL
!        CALL GETCONTROLFPQQ(CONTROL)
!        CONTROL = CONTROL .OR. FPCW$ZERODIVIDE
!        CONTROL = CONTROL .OR. FPCW$INVALID
!        CONTROL = CONTROL .OR. FPCW$OVERFLOW
!        CALL SETCONTROLFPQQ(CONTROL)
!#endif
!$$

! ERROR CHECKING FOR INPUT
IF (NRHS .NE. 1) THEN
    CALL MEXERRMSGTXT('MultMexError: 1 INPUT ARGUMENT IS REQUIRED')
ENDIF

IF (NLHS .NE. 1) THEN
    CALL MEXERRMSGTXT('MultMexError: 1 OUTPUT ARGUMENT IS REQUIRED')
ENDIF

! ASSIGN POINTERS TO THE VARIOUS PARAMETERS
! Input matrix
matAR   =>MXGETPR(PRHS(1))
mA      = MXGETM(PRHS(1))
nA      = MXGETN(PRHS(1))
! Do something meaningful...
classId = mxClassIDFromClassName('double')
plhs(1) = mxCreateDoubleMatrix(mA,nA,0)
AInvR   =>mxGetPr(plhs(1))
call invTest(matAR,mA,nA,AInvR)
! Send the result to the return argument
!(For an alternative way of sending the results to the return arguments - see referenceF90.f90)
END SUBROUTINE MEXFUNCTION

The test function where I do the matrix inversion is as follows:

subroutine invTest(matAR,mA,nA,AInvR)
use mexf90
implicit none

interface
    function dinv(A) result(Ainv)
        use mexf90
        real*8, dimension(:,:), intent(in)      :: A
        real*8, dimension(size(A,1),size(A,2))  :: Ainv
        end function dinv
end interface

integer*4, intent(in)   :: mA, nA
real*8, intent(in)      :: matAR(mA,nA)
real*8, intent(out)     :: AInvR(mA,nA)

AInvR       = dinv(matAR)
end subroutine invTest

and the subroutine which I found on the internet to do the matrix inversion is:

function dinv(A) result(Ainv)
    use mexf90
    real*8, dimension(:,:), intent(in)      :: A
    real*8, dimension(size(A,1),size(A,2))  :: Ainv

    real*8, dimension(size(A,1)) :: work  ! work array for LAPACK
    integer, dimension(size(A,1)) :: ipiv   ! pivot indices
    integer :: n, info, ddd

    ! External procedures defined in LAPACK
    external DGETRF
    external DGETRI

    ! Store A in Ainv to prevent it from being overwritten by LAPACK
    Ainv = A
    n = size(A,1)

    ! DGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call DGETRF(n, n, Ainv, n, ipiv, info)
    if (info /= 0) then
     CALL MEXERRMSGTXT('Matrix is numerically singular!')
    end if

    ! DGETRI computes the inverse of a matrix using the LU factorization
    ! computed by DGETRF.
    call DGETRI(n, Ainv, n, ipiv, work, n, info)

    if (info /= 0) then
     CALL MEXERRMSGTXT('Matrix inversion failed!')
    end if
end function dinv

Answer 1

You've got a mismatch of types for nA and mA in your call to invTest. They should be integer*4 but you have them declared as integer*8 in MEXFUNCTION.

If they need to be integer*8 for the call to mxCreateDoubleMatrix then just make integer*4 copies or possibly try

call invTest(matAR,int(mA),int(nA),AInvR)