Calling LAPACK in MEX files

Hi,

I am trying to write mex files with LAPACK calls, but I keep getting mysterious runtime issues. (MATLAB compiles the mex file, but then while running just crashes with no message.) I looked over my code compared against examples of dgemm and dgesv and I think I'm doing everything right (but there is no dsyev to compare against as far as I'm aware.)

Here is my mex code, which is in projections.c

 #if !defined(_WIN32)
 #define dsyev dsyev_
 #define dgemm dgemm_
 #endif

 #include <stdlib.h>
 #include<stdio.h>
 #include "mex.h"

 #include "lapack.h"
 #include "blas.h"

 void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[]);

 void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
 {

    char *vectors = "N";
    char *uplo = "U";

    double *inMatrix;               
    double *outMatrix; 
    double *eigenvalues;

    mwSize m,n;                 
    int info, lwork;
    int i,j;
    double * work;
    double wkopt;

    inMatrix = mxGetPr(prhs[0]); 
    m = mxGetM(prhs[0]);  
    n = mxGetN(prhs[0]);
    eigenvalues = (double*)malloc(n*sizeof(double));
    plhs[0] = mxCreateDoubleMatrix(m,n,mxREAL);
    outMatrix = mxGetPr(plhs[0]);
    memcpy(outMatrix, inMatrix, m*n*sizeof(double));
    lwork = -1;
    dsyev(vectors, uplo, &n, outMatrix, &n, eigenvalues, &wkopt, &lwork, &info);
  }

and here is my matlab call

 lapacklib = fullfile(matlabroot,'extern','lib',computer('arch'),'microsoft',...
  'libmwlapack.lib');

 mex('-largeArrayDims', '../c/projections.c',lapacklib) %so far so good

 n = 10
 X = randn(n); X = X + X'; 
 Y = projections(X) %here MATLAB crashes

Any idea what I'm doing wrong? I checked to make sure I'm not running an in place operation on the input pointer, but since the matrixdivide.m example does the inplace on the output pointer, I figured that would be ok. Any tips would be appreciated!

Answer 1

[James Tursa]

Some observations:

1) You are using mwSize for the type of integer to use for the LAPACK arguments. Don't do this since this can result in a mismatch between what mwSize is and what the LAPACK/BLAS library functions expect. You should use mwSignedIndex or ptrdiff_t instead. E.g.,

    ptrdiff_t m, n;                 
    ptrdiff_t info, lwork;

2) As written, you have a memory leak for eigenvalues since you use malloc instead of mxMalloc.

3) You are using lwork = -1, which means only a "workspace query". So no eigenvalues will be calculated, only the size of an optimal work array. I assume this was your intent and the only thing you are trying to get working right now is getting that optimal value calculated.

4) I guess I am a little surprised that you don't have to link with libmwblas.lib also. What version of MATLAB are you using? 32-bit or 64-bit?

5) Why do you have a prototype for mexFunction in your code?

6) You use the memcpy function but do not include its header to get the prototype. You should have this in your code:

#include <string.h>

I assume what you have posted is only for the purpose of getting the optimal lwork value, and that you intend to add more code (i.e. another dsyev call) to calculate the eigenvectors/eigenvalues themselves and return them to the workspace.

[ys]

Hi, Observation 1 did turn out to be the bug! Sorry I had resolved it and not replied to the post.

2,3) I was just trying to get it to run without runtime error. Later I did free the malloc and add a second call, so yeah, it's bad code as it stands.

4) I didn't realize I needed to link blas if I'm just using a lapack call? Later I added in blas calls so now it is linked. I'm using 64 bit windows, and it actually worked fine without the blas lib.

5) Didn't realize I didn't need one! Removed!

6) Hmm, thanks, that is a useful tip! I'm trying to keep the code c-friendly in case I eventually port to an all-c-version, which is why I'm avoiding mx functions where possible. But there are definitely some key differences between mex and c!

Thanks for your help and tips!