Showing results for 
Search instead for 
Did you mean: 

creating a MATLAB .mex file using dspsv

I'm trying to use the MKL libraries to build a MATLAB routine to solve a matrix equation with a packed storage symmetric matrix (to save memory). I'm able to get it up and running fine using the following code, but for larger problems (say 5000x5000, which is still relatively small) it will run fine ~50% of the time, and crash MATLAB ~50% of the time. Does anybody see anything that I am obviously doing wrong?
typedef size_t INT;
#define MKL_INT INT
#include "mkl_lapacke.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/* Locals */
lapack_int n, nr, lda, ldb, info;
double *ar, *br;
ar = mxGetPr(prhs[0]); /* pointer to first input matrix */
br = mxGetPr(prhs[1]); /* pointer to second input matrix */
// dimensions of input matrix A //
lda = mxGetM(prhs[0]);
n = mxGetM(prhs[1]);
// dimensions of input matrix B //
ldb = n;
nr = mxGetN(prhs[1]);
if (lda != (ldb*(ldb+1)/2))
mexErrMsgTxt("Error, size mismatch.");
// Local arrays //
lapack_int *ipiv;
ipiv = (lapack_int *) mxMalloc(2*n*sizeof(lapack_int));
//lapack_int ipiv;
info = LAPACKE_dspsv( LAPACK_COL_MAJOR, 'U', n, nr, ar, ipiv, br, ldb);
0 Kudos
2 Replies

The problem appears to be something I'm doing wront in LAPACKE_dspsv. When the only thing I do is remove the size checking for the symmetric matrices and change to dgesv using:
[cpp]info = LAPACKE_dgesv( LAPACK_COL_MAJOR, n, nr, ar, lda, ipiv, br, ldb);[/cpp]
everything seems to work perfectly. At least I haven't been able to crash it yet. Using dspsv, it would generally run once, and then crash the next time it ran. The input matrix is a column matrix that has the right number of elements for a packed storage matrix (n*(n+1)/2), so I have no idea what the problem could be, especially since it seems to work correctly for very small arrays and when it does run it runs fine for large arrays as well. Any ideas?
0 Kudos

Nevermind, I figured it out I think. I was linking to the dynamic libraries and they weren't being unloaded after I ran "clear all". I statically linked and everything seems to be working.
0 Kudos