Get trouble to call Lapack dgesvd with mex

[Jing Qian]

I'm trying to call Lapack subroutine dgesvd from my Fortran code through mex. But no matter what I tried, the dgesvd always return the INFO value as -13. That' means the work space is not correct. I am using the Matab 2010b and Intel Fortran 2013 on Windows 7 64-bit. The libmwlapack library is in Matlab folder C:\Program Files\MATLAB\R2012b\extern\lib\win64\microsoft. Here is my code:

#INCLUDE "FINTRF.H"

SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
MWPOINTER PLHS(*), PRHS(*)
INTEGER NLHS, NRHS, INFO, LWORK
MWPOINTER MXCREATEDOUBLEMATRIX
MWPOINTER MXGETPR
MWPOINTER APR, MPR, NPR, BPR, CPR,WPR
mwSize MX, NX, NN, MM
MWSIGNEDINDEX MXGETM, MXGETN
real*8 :: Mtmp, Ntmp
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: A, C
REAL*8, ALLOCATABLE, DIMENSION(:) :: B
REAL*8, ALLOCATABLE, DIMENSION(:) :: WORK
REAL*8, DIMENSION(1,1) :: VT


APR = MXGETPR(PRHS(1))
MX = mxgetm(prhs(1))
NX = mxgetn(prhs(1))
NN = MIN(MX,NX)
MM = MAX(MX,NX)

ALLOCATE(A(MX,NX))
CALL MXCOPYPTRTOREAL8(APR, A, MX*NX)
ALLOCATE(B(NX))
ALLOCATE(C(MX,NN))


ALLOCATE(WORK(1))
CALL dgesvd('S','N',MX,NX,A,MX,B,C,MX,VT,1,WORK,-1,INFO)
LWORK = INT(WORK(1))
deallocate(WORK)

ALLOCATE(WORK(LWORK))
CALL dgesvd('S','N',MX,NX,A,MX,B,C,MX,VT,1,WORK,LWORK,INFO)

PLHS(1) = MXCREATEDOUBLEMATRIX(MX, NN, 0.0)
CALL MXCOPYREAL8TOPTR(C,PLHS(1),MX*NN)
PLHS(2) = MXCREATEDOUBLEMATRIX(NN, 1, 0.0)
CALL MXCOPYREAL8TOPTR(B,PLHS(2),NN)

RETURN
END SUBROUTINE MEXFUNCTION
Thanks for any help!!!

[Jing Qian]

Actually, I'm using Matlab 2012b. Sorry for mistake.
"Jing Qian" wrote in message <knjfps$fuh$1...@newscl01ah.mathworks.com>...

[James Tursa]

"Jing Qian" wrote in message <knjfps$fuh$1...@newscl01ah.mathworks.com>...

All integers passed to the BLAS & LAPACK routines need to be MWSIGNEDINDEX type. And you need to pay closer attention to the signature of the MATLAB API functions you are calling. E.g.,

> SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
> MWPOINTER PLHS(*), PRHS(*)
> INTEGER NLHS, NRHS, INFO, LWORK

The doc says the signature is:

subroutine mexFunction(nlhs, plhs, nrhs, prhs)
integer*4 nlhs, nrhs
mwPointer plhs(*), prhs(*)

Now, even on a 64-bit system I would typically expect INTEGER to be INTEGER*4, but why depend on this? Also, INFO and LWORK are passed into LAPACK routines, so they definitely should be MWSIGNEDINDEX. E.g., change the above to:

SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
MWPOINTER PLHS(*), PRHS(*)

INTEGER*4 NLHS, NRHS
MWSIGNEDINDEX INFO, LWORK ! probably 8-byte integers on your system

Then consider this declaration:

> MWSIGNEDINDEX MXGETM, MXGETN

The doc says these are supposed to be MWPOINTER (essentially a size_t behind the scenes), so go ahead and use that. E.g.,

MWPOINTER, EXTERNAL :: MXGETM, MXGETN

Then consider this declaration:

> mwSize MX, NX, NN, MM

You use these variables two different ways. One way is as integer arguments to MATLAB API functions like MXCREATEDOUBLEMATRIX (which wants MWSIZE type inputs), and then you also use them in the LAPACK calls to DGESVD (which wants MWSIGNEDINDEX type inputs). So here is another potential for mismatch. Fix this up, probably by using two sets of variables (one typed MWSIZE and another typed MWSIGNEDINDEX).

Then consider this call:

> CALL dgesvd('S','N',MX,NX,A,MX,B,C,MX,VT,1,WORK,-1,INFO)

You pass constants in the arguments (1 and -1). You can get away with this in C/C++ where arguments automatically get promoted to the expected type, but you can't get away with this in Fortran where things are passed by reference. These constant integers are probably 4-bytes each and the routine is probably expecting 8-byte integers. Never use constants as arguments to BLAS or LAPACK routines in Fortran, always use variables. E.g., the above should be something like:

MWSIGNEDINDEX :: PLUS_ONE = 1
MWSIGNEDINDEX :: MINUS_ONE = -1
:
CALL dgesvd('S','N',MX,NX,A,MX,B,C,MX,VT,PLUS_ONE,WORK,MINUS_ONE,INFO)

And for these lines:

> PLHS(1) = MXCREATEDOUBLEMATRIX(MX, NN, 0.0)
> CALL MXCOPYREAL8TOPTR(C,PLHS(1),MX*NN)
> PLHS(2) = MXCREATEDOUBLEMATRIX(NN, 1, 0.0)
> CALL MXCOPYREAL8TOPTR(B,PLHS(2),NN)

just scrub them to replace all constants with variables of the correct type, and make sure the variables you are passing are the correct type as well (e.g., MWSIZE).


James Tursa

[Jing Qian]

Hello, James:
Thank you very much for the great explanation. Following your suggestion, I modified my code and it works now. One more thing I need to point out that the MWSIGNEDINDEX in my machine is corresponding to integer*4. Buy the 64-bit libmwlapack.lib actually accepts integer*8. So, I define the ints used by lapack as integer*8. In addition, the result generated from my code is not same as the matlab function SVD. But manual says SVD does use the function of DGESVD. Do you have any idea about that?

Best regards,
Jing

I post the updated code here in case of someone needs it.

#INCLUDE "FINTRF.H"

SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
MWPOINTER PLHS(*), PRHS(*)
INTEGER NLHS, NRHS

MWPOINTER :: MXCREATEDOUBLEMATRIX
MWPOINTER :: MXGETPR
MWPOINTER, EXTERNAL :: MXGETM, MXGETN
MWPOINTER :: APR, BPR, CPR, WPR
MWSIZE :: M, N, NN, MM
INTEGER*8 :: MX, NX, LWORK, INFO
INTEGER*8 :: ONE = 1
CHARACTER CH1, CH2
DOUBLE PRECISION :: MTMP, NTMP
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: A, C
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: B
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: WORK
DOUBLE PRECISION, DIMENSION(1,1) :: VT

CH1 = 'S'
CH2 = 'N'
M = MXGETM(PRHS(1))
N = MXGETN(PRHS(1))

NN = MIN(MX,NX)
MM = MAX(MX,NX)

MX = M
NX = N

APR = MXGETPR(PRHS(1))
PLHS(1) = MXCREATEDOUBLEMATRIX(M, NN, 0.0)
CPR = MXGETPR(PLHS(1))

PLHS(2) = MXCREATEDOUBLEMATRIX(NN, 1, 0.0)

BPR = MXGETPR(PLHS(2))

LWORK = -1
ALLOCATE(WORK(1))
CALL DGESVD(CH1,CH2,MX,NX,%VAL(APR),MX,%VAL(BPR),
+%VAL(CPR),MX,VT,ONE,WORK,LWORK,INFO)
LWORK = INT(WORK(1))
DEALLOCATE(WORK)

ALLOCATE(WORK(LWORK))
CALL DGESVD(CH1,CH2,MX,NX,%VAL(APR),MX,%VAL(BPR),
+%VAL(CPR),MX,VT,ONE,WORK,LWORK,INFO)
DEALLOCATE(WORK)

RETURN
END

[James Tursa]

"Jing Qian" wrote in message <knlgft$lvk$1...@newscl01ah.mathworks.com>...

> Hello, James:
> Thank you very much for the great explanation. Following your suggestion, I modified my code and it works now. One more thing I need to point out that the MWSIGNEDINDEX in my machine is corresponding to integer*4. Buy the 64-bit libmwlapack.lib actually accepts integer*8. So, I define the ints used by lapack as integer*8.

Hmmm ... that should NEVER happen! In fact, that is the whole point of using MWSIGNEDINDEX so the integer types match up. Can you compile and run this code and post the output:

#INCLUDE "FINTRF.H"
SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
MWPOINTER PLHS(*), PRHS(*)

INTEGER*4 NLHS, NRHS
INTEGER*4, EXTERNAL :: MEXPRINTF
MWSIGNEDINDEX X(2)
MWPOINTER P(2)
MWPOINTER X1, X2
INTEGER*4 K
CHARACTER*50 LINE
X1 = LOC(X(1)) ! you may have to use %LOC(X(1)) here and below, etc
X2 = LOC(X(2))
WRITE(LINE,*) 'Sizeof(MWSIGNEDINDEX) =',X2-X1
K = MEXPRINTF(LINE//ACHAR(10))
X1 = LOC(P(1))
X2 = LOC(P(2))
WRITE(LINE,*) 'Sizeof(MWPOINTER) =',X2-X1
K = MEXPRINTF(LINE//ACHAR(10))
RETURN
END

If you are in fact getting MWPOINTER being 8 bytes and MWSIGNEDINDEX being 4 bytes on a 64-bit machine where the BLAS and LAPACK libraries expect 8 byte integers, I would consider that an error of some type in either your setup, your compiling method, or in the fintrf.h file itself. I will try to look at this on a 64-bit machine I have access to.

> In addition, the result generated from my code is not same as the matlab function SVD. But manual says SVD does use the function of DGESVD. Do you have any idea about that?

No. I am not familiar with the details of how MATLAB calls the LAPACK functions for computing SVD.

> APR = MXGETPR(PRHS(1))
> PLHS(1) = MXCREATEDOUBLEMATRIX(M, NN, 0.0)
> CPR = MXGETPR(PLHS(1))
> PLHS(2) = MXCREATEDOUBLEMATRIX(NN, 1, 0.0)
> BPR = MXGETPR(PLHS(2))

Note, you are still passing constants in the argument lists above, which should be avoided.


James Tursa

[James Tursa]

"James Tursa" wrote in message <knlnf6$hvu$1...@newscl01ah.mathworks.com>...

Can you add these lines to the file as well:

#if defined(MX_COMPAT_32)
K = MEXPRINTF('MX_COMPAT_32 is defined'//ACHAR(10))
#else
K = MEXPRINTF('MX_COMPAT_32 is not defined'//ACHAR(10))
#endif

James Tursa

[Jing Qian]

I got the results like this:

Sizeof(MWSIGNEDINDEX) = 4
Sizeof(MWPOINTER) = 8
MX_COMPAT_32 is defined

what does the MX_COMPAT mean?

"James Tursa" wrote in message <knlos1$mr6$1...@newscl01ah.mathworks.com>...

[James Tursa]

"Jing Qian" wrote in message <knmls5$i0c$1...@newscl01ah.mathworks.com>...

> I got the results like this:
>
> Sizeof(MWSIGNEDINDEX) = 4
> Sizeof(MWPOINTER) = 8
> MX_COMPAT_32 is defined

MX_COMPAT_32 is a macro that, when defined, forces several of the integer types in the MATLAB API (mwSize, mwSignedIndex, etc) to be 4 bytes. I don't know the details of when & how this gets defined in the mex process. On older versions of 32-bit MATLAB I think it is defined by default, and you have to use the -largeArrayDims option in the mex command to get it undefined (which allows 8-byte integers). For later versions of 32-bit MATLAB I think MX_COMPAT_32 is *not* defined by default, meaning that the default is -largeArrayDims and you have to do something to get MX_COMPAT_32 defined. I *thought* that this was true for 64-bit systems as well, i.e. I thought that MX_COMPAT_32 was *not* defined by default on 64-bit systems. But in your setup, for some reason, you have MX_COMPAT_32 defined which is why the BLAS and LAPACK calls are all screwed up when using mwSignedIndex for the integer
types. I will state again that this should not happen regardless of the settings you use. For years I have always thought that mwSignedIndex was set up by TMW to always work for the BLAS & LAPACK interface regardless of how you compiled the code. I now believe this is incorrect, and in fact I also believe the official documented examples that TMW supplies with MATLAB are misleading or wrong. While there may be separate 4-byte integer and 8-byte integer versions of API functions such as mxCreateDoubleMatrix etc, I don't think that is true of the BLAS and LAPACK libraries (I may have to research this).

Bottom line for me is that in the future I will always use mwPointer (in Fortran) for the integer types in BLAS and LAPACK routine calls. Thanks for your post and for running the test for me ... I learned something new.

Getting back to why you have MX_COMPAT_32 defined on your 64-bit system, I will ask you exactly what command(s) you are using to compile the mex routine? What are you doing that is causing the macro MX_COMPAT_32 to be defined?

James Tursa

[James Tursa]

"James Tursa" wrote in message <knn1bu$j9p$1...@newscl01ah.mathworks.com>...
:

> While there may be separate 4-byte integer and 8-byte integer versions of API functions such as mxCreateDoubleMatrix etc, I don't think that is true of the BLAS and LAPACK libraries (I may have to research this).

Confirmed. I found this statement in the doc:

"Use the -largeArrayDims option; the mwlapack and mwblas libraries only support 64-bit integers for matrix dimensions."

So mwSignedIndex should *not* be used for BLAS and LAPACK calls. You should always use mwPointer.

James Tursa

[James Tursa]

"James Tursa" wrote in message <knn1bu$j9p$1...@newscl01ah.mathworks.com>...

>
> Bottom line for me is that in the future I will always use mwPointer (in Fortran) for the integer types in BLAS and LAPACK routine calls. Thanks for your post and for running the test for me ... I learned something new.

After thinking about it some more I realize that isn't always going to work either. I will need to research when the BLAS and LAPACK libraries switched to 8-byte integer types.

James Tursa