Using LAPACK and BLAS Functions :: Creating C Language MEX-Files (External Interfaces)

External Interfaces

Using LAPACK and BLAS Functions

LAPACK is a large, multiauthor Fortran subroutine library that MATLAB uses for numerical linear algebra. BLAS, which stands for Basic Linear Algebra Subroutines, is used by MATLAB to speed up matrix multiplication and the LAPACK routines themselves. The functions provided by LAPACK and BLAS can also be called directly from within your C MEX-files.

This section explains how to write and build MEX-files that call LAPACK and BLAS functions. It provides information on

Specifying the Function Name

When calling an LAPACK or BLAS function, some platforms require an underscore character following the function name in the call statement.

On the PC and HP platforms, use the function name alone, with no trailing underscore. For example, to call the LAPACK dgemm function, use

```
dgemm (arg1, arg2, ..., argn);
```

On the LINUX, Solaris, and Macintosh platforms, add the underscore after the function name. For example, to call dgemm on any of these platforms, use

```
dgemm_ (arg1, arg2, ..., argn);
```

Calling LAPACK and BLAS Functions from C

Since the LAPACK and BLAS functions are written in Fortran, arguments passed to and from these functions must be passed by reference. The following example calls dgemm, passing all arguments by reference. An ampersand (&) precedes each argument unless that argument is already a reference.

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, mxArray 
*prhs[])
{
  double *A, *B, *C, one = 1.0, zero = 0.0;
  int m,n,p; 
  char *chn = "N";

  A = mxGetPr(prhs[0]);
  B = mxGetPr(prhs[1]);
  m = mxGetM(prhs[0]);
  p = mxGetN(prhs[0]);
  n = mxGetN(prhs[1]);

  if (p != mxGetM(prhs[1])) {
    mexErrMsgTxt("Inner dimensions of matrix multiply do not
        match");
  }

  plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
  C = mxGetPr(plhs[0]);

  /* Pass all arguments to Fortran by reference */
  dgemm (chn, chn, &m, &n, &p, &one, A, &m, B, &p, &zero, C, &m);
}

Handling Complex Numbers

MATLAB stores complex numbers differently than Fortran. MATLAB stores the real and imaginary parts of a complex number in separate, equal length vectors, pr and pi. Fortran stores the same number in one location with the real and imaginary parts interleaved.

As a result, complex variables exchanged between MATLAB and the Fortran functions in LAPACK and BLAS are incompatible. MATLAB provides conversion routines that change the storage format of complex numbers to address this incompatibility.

Input Arguments. For all complex variables passed as input arguments to a Fortran function, you need to convert the storage of the MATLAB variable to be compatible with the Fortran function. Use the mat2fort function for this. See the example that follows.

Output Arguments. For all complex variables passed as output arguments to a Fortran function, you need to do the following:

When allocating storage for the complex variable, allocate a real variable with twice as much space as you would for a MATLAB variable of the same size. You need to do this because the returned variable uses the Fortran format, which takes twice the space. See the allocation of zout in the example that follows.
Once the variable is returned to MATLAB, convert its storage so that it is compatible with MATLAB. Use the fort2mat function for this.

Example - Passing Complex Variables. The example below shows how to call an LAPACK function from MATLAB, passing complex prhs[0] as input and receiving complex plhs[0] as output. Temporary variables zin and zout are used to hold prhs[0] and plhs[0] in Fortran format.

#include "mex.h" 
#include "fort.h"      /* defines mat2fort and fort2mat */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, mxArray 
*prhs[])
{ 
  int lda, n; 
  double *zin, *zout;
  lda = mxGetM(prhs[0]); 
  n = mxGetN(prhs[0]); 

  /* Convert input to Fortran format */
  zin = mat2fort(prhs[0], lda, n);

  /* Allocate a real, complex, lda-by-n variable to store output 
*/
  zout = mxCalloc(2*lda*n, sizeof(double));

  /* Call complex LAPACK function */ 
  zlapack_function(zin, &lda, &n, zout);

  /* Convert output to MATLAB format */ 
  plhs[0] = fort2mat(zout, lda, lda, n);

  /* Free intermediate Fortran format arrays */ 
  mxFree(zin);
  mxFree(zout); 
}

Preserving Input Values from Modification

Many LAPACK and BLAS functions modify the values of arguments passed in to them. It is advisable to make a copy of arguments that can be modified prior to passing them to the function. For complex inputs, this point is moot since the mat2fort version of the input is a new piece of memory, but for real data this is not the case.

The following example calls an LAPACK function that modifies the first input argument. The code in this example makes a copy of prhs[0], and then passes the copy to the LAPACK function to preserve the contents of prhs[0].

/* lapack_function modifies A so make a copy of the input */ 
m = mxGetM(prhs[0]); 
n = mxGetN(prhs[0]); 
A = mxCalloc(m*n, sizeof(double)); 

/* Copy mxGetPr(prhs[0]) into A */ 
temp = mxGetPr(prhs[0]); 
for (k = 0; k < m*n; k++) { 
   A[k] = temp[k]; 
   }

/* lapack_function does not modify B so it is OK to use the input 
directly */ 
B = mxGetPr(prhs[1]);
lapack_function(A, B);       /* modifies A but not B */

/* Free A when you are done with it */ 
mxFree(A);

Building the C MEX-File

The examples in this section show how to compile and link a C MEX file, myCmexFile.c, on the platforms supported by MATLAB. In each example, the term <matlab> stands for the MATLAB root directory.

Building on the PC. If you build your C MEX-file on a PC platform, you need to explicitly specify a library file to link with.

On the PC, use this command if you are using the Lcc compiler that ships with MATLAB:

mex myCmexFile.c <matlab>/extern/lib/win32/lcc/libmwlapack.lib

Or, use this command if you are using Microsoft Visual C++ as your C compiler:

mex myCmexFile.c 
<matlab>/extern/lib/win32/microsoft/msvc60/libmwlapack.lib

Building on Other Platforms. On all other platforms, you can build your MEX-file as you would any other C MEX-file. For example,

```
mex myCmexFile.c
```

MEX-Files That Perform Complex Number Conversion. MATLAB supplies the files fort.c and fort.h, which provide routines for conversion between MATLAB and FORTRAN complex data structures. These files define the mat2fort and fort2mat routines mentioned previously under Handling Complex Numbers.

If your program uses these routines, then you need to:

Include the fort.h file in your program, using, #include "fort.h". See the example above.
Build the fort.c file with your program. Specify the pathname, <matlab>/extern/examples/refbook for both fort.c and fort.h in the build command, (where <matlab> stands for the MATLAB root directory).

On the PC, use either one of the following.

mex myCmexFile.c <matlab>/extern/examples/refbook/fort.c 
-I<matlab>/extern/examples/refbook 
<matlab>/extern/lib/win32/microsoft/msvc60/libmwlapack.lib
mex myCmexFile.c <matlab>/extern/examples/refbook/fort.c 
-I<matlab>/extern/examples/refbook 
<matlab>/extern/lib/win32/lcc/libmwlapack.lib

For all other platforms, use

mex myCmexFile.c <matlab>/extern/examples/refbook/fort.c 
-I<matlab>/extern/examples/refbook

Example - Symmetric Indefinite Factorization Using LAPACK

You will find an example C MEX-file that calls two LAPACK functions in the directory <matlab>/extern/examples/refbook, where <matlab> stands for the MATLAB root directory. There are two versions of this file:

utdu_slv.c - calls functions zhesvx and dsysvx, and thus is compatible with the PC and HP platforms.
utdu_slv_.c - calls functions zhesvx_ and dsysvx_, and thus is compatible with the LINUX and Solaris platforms.

Calling LAPACK and BLAS Functions from Fortran

You can make calls to the LAPACK and BLAS functions used by MATLAB from your Fortran MEX files. The following is an example program that takes two matrices and multiplies them by calling the LAPACK routine, dgemm:

      subroutine mexFunction(nlhs, plhs, nrhs, prhs)
      integer plhs(*), prhs(*)
      integer nlhs, nrhs
      integer mxcreatedoublematrix, mxgetpr
      integer mxgetm, mxgetn
      integer m, n, p
      integer A, B, C
      double precision one, zero, ar, br
      character ch1, ch2

      ch1 = 'N'
      ch2 = 'N'
      one = 1.0
      zero = 0.0

      A = mxgetpr(prhs(1))
      B = mxgetpr(prhs(2))
      m = mxgetm(prhs(1))
      p = mxgetn(prhs(1))
      n = mxgetn(prhs(2))

      plhs(1) = mxcreatedoublematrix(m, n, 0.0)
      C = mxgetpr(plhs(1))
      call mxcopyptrtoreal8(A, ar, 1)
      call mxcopyptrtoreal8(B, br, 1)

      call dgemm (ch1, ch2, m, n, p, one, %val(A), m,
     +            %val(B), p, zero, %val(C), m)

      return
      end

Building the Fortran MEX-File

The examples in this section show how to compile and link a Fortran MEX file, myFortranmexFile.F, on the platforms supported by MATLAB. In each example, the term <matlab> stands for the MATLAB root directory.

Building on the PC. On the PC, using Visual Fortran, you will have to link against a library called libdflapack.lib:

mex -v myFortranMexFile.F 
<matlab>/extern/lib/win32/digital/df60/libdflapack.lib

Building on Other UNIX Platforms. On the UNIX platforms, you create the MEX file as follows:

```
mex -v myFortranMexFile.F
```

Large File I/O Debugging C Language MEX-Files