Matlab crashed when running Pardiso mex file with certain input size

Hi everybody

I am experienced in Matlab, but not that much in C. I work on this issue for a while now and I hope someone can point out the problem.

I want to use PARDISO to solve a large system of equations (~5e5 x 5e5 sparse matrix with ~5e6 non-zeros). Currently I apply the Matlab's backslash solver. I wrote a mex file to call pardiso and was successful in compliling it with the Visual Studio C library. I tested several smaller linear systems (8x8) with type=11 (unsymmetric, real) and they worked fine. However, Matlab sooner or later crashes, either by repeating the mex function, or by increasing the matrix size. I am pretty sure there is a memory problem, but I cannot locate it. Also, I don't understand any of the temporary crash log files. I tried around with int32 and int64, modifying the input integers in both Matlab and C, but the default option is the only one that works for small matrices.

My computer has 128 GB memory and it does not show any specific peak of memory usage when it crashes.

I compiled the c-file with the following command (it compiled correctly, and I could run it for small matrices, so I think that is ok):

mex pardiso_mex_corr.c -I"C:\Program Files (x86)\Intel\oneAPI\mkl\2025.1\include" "C:\Program Files (x86)\Intel\oneAPI\mkl\2025.1\lib\mkl_rt.lib"

I defined the input of the mex function as arrays in CSR format with double (A,B) and as int32 (JA,IA,na) in Matlab. Since the function gives correct results, I suppose that thisis ok.

Here is what I did:

#include "mex.h"
#include "mkl.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    /* Check number of inputs and outputs */
    if (nrhs != 5) {
        mexErrMsgIdAndTxt("MATLAB:PARDISO:invalidNumInputs", "Five input arguments required.");
    }

    /* Get input arrays from MATLAB */
    double *A = (double*)mxGetPr(prhs[0]);   // Matrix A
    int *JA = (int*)mxGetData(prhs[1]);  // Column indices
    int *IA = (int*)mxGetData(prhs[2]);  // Row indices
    double *B = mxGetPr(prhs[3]);   // RHS vector B
    int na = (int)mxGetScalar(prhs[4]);   // Matrix size
    
    plhs[0] = mxCreateDoubleMatrix(na, 1, mxREAL);  // Create a MATLAB double column vector
    double *xx = mxGetPr(plhs[0]);                // Pointer to the data (to be filled by PARDISO)    

    /* Setup PARDISO parameters */
    int i;
    int mtype = 11; /* Real unsymmetric matrix */
    void *pt[64] = {NULL};   /* Internal solver memory pointer */
    int iparm[64];
    int maxfct = 1, mnum = 1, phase, error, msglvl = 1;
    double ddum;    // Double dummy
    int idum;       // Integer dummy

    /* Initialize PARDISO control parameters */
    for (i = 0; i < 64; i++) {
        iparm[i] = 0;
    }
    iparm[0] = 1;  /* No solver default */
    iparm[1] = 0;  /* If set to 2: Fill-in reordering from METIS */
    iparm[2] = 1;  /* Number of processors, value of OMP_NUM_THREADS */
    iparm[3] = 0;  /* No iterative-direct algorithm */
    iparm[7] = 20; /* Max numbers of iterative refinement steps */
    iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
    iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
    iparm[18] = -1; /* Output: Mflops for LU factorization */
    iparm[19] = 0;  /* Output: Numbers of CG Iterations */

    /* Symbolic factorization */
    phase = 11;
    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &na, A, IA, JA,
            &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
    if (error != 0) {
        mexErrMsgIdAndTxt("MATLAB:PARDISO:symbolicFactorizationError", "Error during symbolic factorization: %d", error);
    }

    /* Numerical factorization */
    phase = 22;
    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &na, A, IA, JA,
            &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
    if (error != 0) {
        mexErrMsgIdAndTxt("MATLAB:PARDISO:numericalFactorizationError", "Error during numerical factorization: %d", error);
    }

    /* Back substitution and iterative refinement */
    phase = 33;
    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &na, A, IA, JA,
            &idum, &nrhs, iparm, &msglvl, B, xx, &error);
    if (error != 0) {
        mexErrMsgIdAndTxt("MATLAB:PARDISO:backSubstitutionError", "Error during back substitution: %d", error);
    }

    if (nrhs != 5) {
        mexErrMsgIdAndTxt("MATLAB:PARDISO:invalidNumInputs", "Five input arguments required.");
    }

    /* Terminate and release solver memory */
    phase = -1;
    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &na, &ddum, IA, JA,
            &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

}

[JONAS BRUNO]

One thing I observed is that for small (8x8) sparse systems that are symmetric, it works fine, but for non symmetric ones, results start to deviate from the accurate solution and after running it 4, 5 times, matlab crashes.

Answer 1

[James Tursa]

I haven't had time to look at this thoroughly, and I am unfamiliar with PARDISO, but the symptoms you describe are usually the case of mismatched argument types or incorrect/incomplete argument allocations, which both lead to accessing invalid memory downstream and MATLAB crashes. Particularly the fact that small problems seem to work fine (the invalid memory wasn't enough to be critical to running) but large problems don't. My advice is to double check that each and every argument variable is the correct type and arrays are allocated to the proper size. E.g.,

It appears to me that nrhs should be an input to PARDISO? I.e., the "number of right hand sides to be solved for". But that conflicts with the use of nrhs in the mexFunction argument list. Two different things. So it appears to me you have a mismatch here, and nrhs for the PARDISO call needs to be something different. E.g., maybe a different named variable such as pnrhs = 1; (I kinda suspect this is the root cause of your crash)

I don't know how to reconcile the signature of the first argument to PARDISO. The doc states both this:

_MKL_DSS_HANDLE_t pt

and this:

pt

Array with size of 64.

That is, typically when I see the first part above in the function signature I would put in a line of code to exactly match that to declare pt:

_MKL_DSS_HANDLE_t pt;

and then pass pt as the first parameter in the call. So what does the "Array with size of 64" comment mean later on in the doc? Is pt supposed to be declared as an array of _MKL_DSS_HANDLE_t like this:

_MKL_DSS_HANDLE_t pt[64] = {0};

Or ...???

This is confusing to me. Maybe looking at what _MKL_DSS_HANDLE_t actually is behind the scenes will clear this up.

Are the arrays really the size that PARDISO is expecting? E.g., the doc calls out some arguments as being arrays of size n*nrhs. Double check that all the arrays in the arguments are the correct size and type. Etc.

Are the integer sizes really what PARDISO is expecting? I.e., is int really the same size as MLK_INT? You should check that sizeof(int) matches sizeof(MLK_INT) and mxGetElementSize(prhs[1]) and mxGetElementSize(prhs[2]) on entry to the mex function. And that the inputs are integer arrays. Frankly, in addition to this check, I would actually use MLK_INT as the explicit type for all integer argument variables to PARDISO.

Bottom line, put in every single input argument check you can think of up front to make your mex function robust. Is prhs[0] really a full real double matrix? Is prhs[3] really a full real double vector of the correct size? Is prhs[4] really numeric and has a reasonable value? Etc.

[JONAS BRUNO]

Dear James

Many thanks for your comments! I will go through all you suggestions. I hope that solves the problem.

Thanks again

J.