MATLAB Answers

Inverse matrix in fortran mexFunction - matlab lapack library

37 views (last 30 days)
Tue
Tue on 21 Jun 2013
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

  1 Comment

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

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 24 Jun 2013
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.

  1 Comment

Sign in to comment.

More Answers (1)

zohar
zohar on 24 Jun 2013
Edited: zohar on 24 Jun 2013
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 !?!?

  2 Comments

Tue
Tue on 24 Jun 2013
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
Muthu Annamalai on 24 Jun 2013
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.

Sign in to comment.

Sign in to answer this question.