using BLAS function in a computational routine in mex

I am using Matlab R2013a 64bit. I am writing a mex function in fortran. I am using the Intel Visual Fortran 14.0 compiler. I am calling a BLAS function. I have the following problem: when I call the BLAS function from the main body of the mex function, it works; but when I call it from a different routine, it crashes Matlab. The following is a minimal code showing this problem, calling the BLAS function dnrm2.

#include "fintrf.h"

!============================================================================

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

implicit none

mwPointer plhs(*), prhs(*)

integer nlhs, nrhs

mwPointer mxCreateDoubleScalar, mxGetPr, mxGetM

integer n

mwPointer x

double precision norm1x, norm2x

double precision dnrm2, compute_norm

integer complexFlag

complexFlag = 0

n = int4(mxGetM(prhs(1)))

x = mxGetPr(prhs(1))

norm1x = dnrm2(n, %val(x), 1) ! WORKS

plhs(1) = mxCreateDoubleScalar(norm1x)

norm2x = compute_norm(n, %val(x)) ! DOES NOT WORK

plhs(2) = mxCreateDoubleScalar(norm2x)

end subroutine mexFunction

!============================================================================

function compute_norm(n, x)

implicit none

integer n

double precision x(n)

double precision compute_norm, dnrm2

compute_norm = dnrm2(n, x, 1)

end function compute_norm

To build, I use

mex -v -largeArrayDims mynorm.F90 'C:\Program Files\MATLAB\R2013a\extern\lib\win64\microsoft\libmwblas.lib'

I am wondering what might be the reason. Thank you.

Answer 1

[James Tursa]

For -largeArrayDims I might expect the BLAS integer arguments to be 64-bit. You are passing 32-bit constant 1's in the 3rd argument. You can get away with this in C/C++ where arguments are automatically type promoted, but you can't get away with this in Fortran if there is a mismatch since there is no automatic type promotion in Fortran. In general, always use variables to pass values to BLAS routines so your code can be robust across different platforms. I would try this as a first guess at a fix:

mwSignedIndex :: ONE = 1
    :
norm1x = dnrm2(n, %val(x), ONE)

and similarly in the compute_norm function:

mwSignedIndex :: ONE = 1
   :
compute_norm = dnrm2(n, x, ONE)

As to why the dnrm2 call in mexFunction worked whereas the dnrm2 call in compute_norm crashed, it may be just luck based on the current state of stack memory because of the calling sequence just prior to the dnrm2 call.

Another possibility is that you only called the mex routine with 1 requested output. In that case the plhs array is only sized for 1 return value, and plhs(2) is writing off the end of the array which will seg fault. You should have protections in your code for this. E.g.,

if( nlhs >= 2 )
    norm2x = compute_norm(n, %val(x)) ! DOES NOT WORK
    plhs(2) = mxCreateDoubleScalar(norm2x)
endif

You should also have something at the beginning of your code to check the input argument. E.g.,

if( nrhs /= 1 ) then
    call mexErrMsgTxt("Need exactly one double vector input")
endif
if( mxIsDouble(prhs(1)) == 0 ) then
    call mexErrMsgTxt("Need exactly one double vector input")
endif
if( mxGetM(prhs(1)) /= 1 .AND. mxGetN(prhs(1)) /= 1 ) then
    call mexErrMsgTxt("Need exactly one double vector input")
endif

[Siva]

Doing

mwSignedIndex :: one = 1  
compute_norm = dnrm2(n, x, one)

solved the problem. Thank you.

I have a follow up question.

I also changed the declaration of n from

integer n

to

mwIndex n

and that works also.

I am wondering the following. 1. How could it work for both declarations of n? 2. I did

>> version('-blas')

and I got

Intel® Math Kernel Library Version 10.3.11 Product Build 20120606 for Intel® 64 architecture applications

The Intel MKL manual lists the types of both the arguments n and incx as INTEGER (I think independently of the platform). But mwIndex and mwSignedIndex are both INTEGER*8. How does this work? In fact, declaring

integer :: one = 1

does not work.

[James Tursa]

You need to match the signature of the particular BLAS you are using. For MATLAB supplied BLAS libraries, that is what mwSignedIndex is for. If you are linking to a different BLAS, then mwSignedIndex may not be appropriate and you may need to use something else as specified in that BLAS doc.

Note that "integer" in Fortran is the default integer, which could be 32-bit or 64-bit. Although not guaranteed, my experience has been that integer is typically a 32-bit quantity on both 32-bit and 64-bit installations, unless special compiler flags are used.

As far as the integer vs mwIndex is concerned, the answer depends. Passing the address of a 32-bit integer to a routine expecting a 64-bit integer will often result in a seg fault because of memory access violations. But passing the address of a 64-bit integer to a routine expecting a 32-bit integer can sometimes work depending on the byte ordering of the integers on the particular machine involved, and whether the value actually fits in the lower 32-bits.

Regardless, the advice is the same ... always use the signature for the routine that is listed in the doc. If you want your function to work for multiple choices of BLAS library where some might use 32-bit integers and some might use 64-bit integers, then you will have to do something special in your code to account for this. You will not be able to blindly link with whatever BLAS you want and expect things to work.

[Siva]

Thank you. Just to clarify again, if I use MATLAB-supplied BLAS, can I always use mwSignedIndex or mwIndex to declare variables that are meant to denote sizes of arrays or indices into arrays without worrying about the signature in the BLAS documentation?

[James Tursa]

If you link with the MATLAB supplied BLAS library, you should always be able to use mwSignedIndex for the integer arguments. What the mwSignedIndex maps to (integer*4 or integer*8) is supposed to all happen automatically for you correctly via the mex command depending on whether your installation is 32-bit or 64-bit and the largeArrayDims flag setting.

That being said, I have not personally tried all of these combinations. There may be some combination that doesn't work that I am not aware of, in which case one would have to manually set the integer argument types directly instead of using the mwSignedIndex macro.

If you link with a 3rd party BLAS library, it is probably best to hard code the integer type depending on that 3rd party doc, and not use mwSignedIndex.

[Siva]

Thank you. That clarifies it.