Trouble Calling lapack functions from mex files on 64 bit Windows

Hello, all,

I am having issues when running the following code, which is supposed to convert the coefficients of a taylor series approximation of order n to those of a pade approximant with numerator and denominator both of order n/2. The file compiles with no errors using the -largeArrayDims flag; however, when I run the file on an example, it produces the appropriate output and then matlab closes. I have also run into issues with other variables in the parent function which calls this routine changing. For instance, I have an array of random indices in the parent function, which is not an input to this mex routine, that gets altered when I make the call to the mex routine. Any help would be very much appreciated. Thanks!

Jonah

#if !defined(WIN32) #define dgesv dgesv #endif

#include "mex.h" #include "lapack.h" #include "matrix.h"

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

    /* initialize variables */
    double *Tlhs, *trhs;
    double *Sr, *Si;
    mxArray *mxPivot;
    mwSignedIndex *iPivot;
    mwSignedIndex dims[2], info;

    double *ar, *ai;
    double *nr, *ni, *dr, *di;
    int ii, jj, nc;
    mwSignedIndex nn;
    mwSignedIndex one=1;

    /* Get input variables */
    ar=mxGetPr(prhs[0]);
    ai=mxGetPr(prhs[1]);

    /* Find the number of Pade coefficients in the num and den */
    nc=(max((int)mxGetM(prhs[0]),(int)mxGetN(prhs[0]))+1)/2;
    dims[0]=(mwSignedIndex)(2*nc-2);
    dims[1]=(mwSignedIndex)(2*nc-2);

    /* Allocate memory to the output array */
    plhs[0]=mxCreateDoubleMatrix(nc,1,mxREAL);
    plhs[1]=mxCreateDoubleMatrix(nc,1,mxREAL);
    plhs[2]=mxCreateDoubleMatrix(nc,1,mxREAL);
    plhs[3]=mxCreateDoubleMatrix(nc,1,mxREAL);

    /* Assign the real and imaginary output ptrs */
    nr=(double *)mxGetPr(plhs[0]);
    ni=(double *)mxGetPr(plhs[1]);
    dr=(double *)mxGetPr(plhs[2]);
    di=(double *)mxGetPr(plhs[3]);
    dr[0]=1;
    di[0]=0;    

    Tlhs = malloc((2*nc-2)*(2*nc-2)*sizeof(ar[0]));
    trhs = malloc((2*nc-2)*sizeof(ar[0]));
    Sr = malloc(nc*sizeof(ar[0]));
    Si = malloc(nc*sizeof(ar[0]));

    /* Construct the Toeplitz Matrix T */
    for(ii=0;ii<nc;ii++)
    {
        for(jj=0;jj<nc;jj++)
        {
            if(jj<nc-1 && ii<nc-1)
            {
                Tlhs[ii+jj*(2*nc-2)]=ar[nc-1-jj+ii];
                Tlhs[nc-1+ii+(nc-1+jj)*(2*nc-2)]=ar[nc-1-jj+ii];
                Tlhs[ii+(nc-1+jj)*(2*nc-2)]=-ai[nc-1-jj+ii];
                Tlhs[nc-1+ii+jj*(2*nc-2)]=ai[nc-1-jj+ii];
            };
        };
        Sr[ii]=ar[ii];
        Si[ii]=ai[ii];
        if(ii<nc-1)
        {
            trhs[ii]=ar[nc+ii];
            trhs[nc-1+ii]=ai[nc+ii];
        };
    };
      
    /* Create inputs for DGESV */
    mxPivot = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
    iPivot = (mwSignedIndex*)mxGetData(mxPivot);

    /* Call LAPACK */
    nn=(2*nc-2);
    dgesv(&nn,&one,Tlhs,&nn,iPivot,trhs,&nn,&info);

    for(ii=0;ii<nc;ii++)
    {
        if(ii<nc-1)
        {
            dr[ii+1]=-trhs[ii];
            di[ii+1]=-trhs[nc-1+ii];
        }
        nr[ii]=0;
        ni[ii]=0;
        for(jj=0;jj<=ii;jj++)
        {
            nr[ii]=nr[ii]+Sr[ii-jj]*dr[jj]-Si[ii-jj]*di[jj];
            ni[ii]=ni[ii]+Sr[ii-jj]*di[jj]+Si[ii-jj]*dr[jj];
        };
    };
  
    free(Sr);
    free(Si);
    free(Tlhs);
    free(trhs);
    mxDestroyArray(mxPivot);

};

Answer 1

[James Tursa]

Why are you using an mxINT32_CLASS for creating mxPivot? Seems like this would only work if mwSignedIndex is a 32-bit integer on your 64-bit system, which it probably isn't. In fact, why create an mxArray at all here, since you are only creating a temporary array for dgesv to use, not a variable to return back to MATLAB? So I would try changing these lines:

    mxPivot = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
    iPivot = (mwSignedIndex*)mxGetData(mxPivot);
        :
    mxDestroyArray(mxPivot);

to these:

    iPivot = (mwSignedIndex*)mxMalloc(dims[0]*dims[1]*sizeof(*iPivot));
        :
    mxFree(iPivot);

[Jonah A Reeger]

James,

Thank you. This seems to have corrected the issue. I had even tried using mxINT64_CLASS in the mxCreateNumericArray call, but never with any luck. It seems that now both of my issues have been taken care of.

Jonah