segfault in matrix inversion using LAPACK in mex

 0 down vote favorite

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

[James Tursa]

Which BLAS and LAPACK are you linking with? For 64-bit MATLAB supplied libraries, I might have expected these routines to require 8-byte integers for the arguments instead of 4-byte integers. E.g., what happens if you try this in the dinv function:

      integer*8, dimension(size(A,1)) :: ipiv   ! pivot indices
      integer*8 :: n, info, ddd

Also, another integer*4 vs integer*8 issue is this:

integer*8               :: err, nA, mA
    :
call invTest(matAR,mA,nA,AInvR)
    :
subroutine invTest(matAR,mA,nA,AInvR)
    :
integer*4, intent(in)   :: mA, nA

So change that last line to:

integer*8, intent(in)   :: mA, nA

[Ugur]

Many thanks, that solved the problem. It was the 32 bit integer declaration in "dinv" function which was causing this error.