MEX: MATLAB crashes when filling array

As my first attempt on a mex file (and c programming in general) I wanted to implement a routine that computes the potential from a spherical harmonic model, at a number of coordinate pairs (lat,lon,h). Basically, this means that I am feeding the mex routine two large matrices (2191x2191), called C and S, along with three vectors (42x1), called lat, lon and h. At each point I then have to loop over all my coefficients in order to compute the potential value. For each coefficient I also have to compute a Legendre value, which I am storing in a similar size (2191x2191) matrix called Pnm. Furthermore, in order to compute the Legendre values, I am computing a vector of seed values (1x2191). Below is an outline of the mex/c file:

#include "mex.h"
#include <math.h>

/* Define constant variables */
const double pi = 3.14159265358979323846;
const double deg2rad = 0.017453292519943;
const double wgs84a = 6378137.0;
const double wgs84e = .081819190842622;
const double sqrt3 = 1.732050807568877;
const double isqrt2 = 0.707106781186547;

/* Function declarations */
double gravpot(double *val, double *lat, double *lon, double *h_ell, int nlat,
               int nlon, int *nmax, int *mmax, double *C, double *S,
               double *shGM, double *sha);

/* ===========================================================================
 * ===== MEX GATEWAY FUNCTION ================================================
 * ======================================================================== */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs,
       const mxArray *prhs[])
{

    /* Pointers to input and output arguments */
    double *xval, *yval, *zval, *lat, *lon, *h_ell, *C, *S, *shGM, *sha;
    unsigned char *lcomp;
    int *nmax, *mmax;

    /* Other */
    int i, j, nlon, nlat;

    /* Check input arguments ------------------------------------------------ */

    /* Extract input data --------------------------------------------------- */

    /* Assign pointers to input */
    lcomp = (unsigned char *) mxGetData(prhs[0]);
    lat   = mxGetPr(prhs[1]);
    lon   = mxGetPr(prhs[2]);
    h_ell = mxGetPr(prhs[3]);
    C     = mxGetPr(prhs[4]);
    S     = mxGetPr(prhs[5]);
    shGM  = mxGetPr(prhs[6]);
    sha   = mxGetPr(prhs[7]);
    nmax  = (int*) mxGetData(prhs[8]);
    mmax  = (int*) mxGetData(prhs[9]);

    /* Get dimensions of input */
    nlat = mxGetM(prhs[1]);
    nlon = mxGetN(prhs[1]);

    /* Change units of lat and lon */
    for (i = 0; i < nlat; i++) {
      for (j = 0; j < nlon; j++) {
        lat[i + j * nlat] = lat[i + j * nlat] * deg2rad;
        lon[i + j * nlat] = lon[i + j * nlat] * deg2rad;
      }
    }

    /* Call computation routine --------------------------------------------- */

    if (*lcomp == 1) {

      /* Create matrix for the return argument */
      plhs[0] = mxCreateDoubleMatrix(nlat,nlon,mxREAL);

      /* Assign pointers to output, this matrix is initialised as zeros */
      xval = mxGetPr(plhs[0]);

      /* Compute gravitational potential */
      gravpot(xval, lat, lon, h_ell, nlat, nlon, nmax, mmax, C, S, shGM, sha);

    }

    /* ---------------------------------------------------------------------- */

}

/* ===========================================================================
 * ===== GRAVITATIONAL POTENTIAL =============================================
 * ======================================================================== */
double gravpot(double *val, double *lat, double *lon, double *h_ell,
               int nlat, int nlon, int *nmax, int *mmax, double *C, double *S,
               double *shGM, double *sha)
{
  /* Declare variables */
  const double sfac = pow(10, -280), ifac = pow(10, 280);
  const int mid = (*nmax+1)*(*nmax+2)-1, mid2 = (*nmax+1)*(*nmax+1)-1;
  double theta, lambda, h, x, z, r, R_N, sphlat;
  double Pnm[mid2+1], Ps[*nmax+1], coef[mid2+1];
  double gnm, hnm, GMr, ar, arn, nsum1, nsum2, msum, t, u, u2;
  double P, Pp1, Pp2, loncos, lonsin;
  int i, j, n, m, id;

    /* Compute sectorial values of associated Legendre polynomials ========== */

    /* Define seed values ( divided by u^m ) */
    Ps[0] = 1.0 * sfac;
    Ps[1] = sqrt3 * sfac;

    /* Compute sectorial values, eq. 13 ( divided by u^m ) */
    for (n = 2; n <= *nmax; n++) { 
      Ps[n] = sqrt( (2.0*n+1) / (2.0*n) ) * Ps[n-1];
    }

    /* ====================================================================== */

    /* NEW SECTION: Filling coef matrix ===================================== */

    coef[0] = 1.0;

    /* ====================================================================== */

    /* Loop through computation points ====================================== */

    /* Loop over latitude */
    for (i = 0; i < nlat; i++) {

      /* Loop over degree */
      for (n = 0; n <= *nmax; n++) {

        /* Loop over order */
        for (m = 0; m <= *mmax; m++) {

          Pnm[n+m*(*nmax+1)] = 1.0;

        }

      }

      /* Loop over longitude */
      for (j = 0; j < nlon; j++) {
    
        val[i+j*nlat] = 1.0;

      }

    }

    /* ====================================================================== */

}

The code above has been stripped from a lot of contents, but the issue withstands for the above piece of code. nnamx = mmax = 2191 are integer numbers. shGM and sha are scalar doubles precision numbers. lcomp is an integer number. If I comment out the line:

coef[0] = 1.0;

There are no problems, but if I attempt to fill the array "coef", MATLAB crashes. I have no problems compiling the file, the problem comes when calling the function. MATLAB throws no error message, but simply crashes. Can anyone tell me if I am doing something wrong? Or am I simply hitting some memory limit?

The code is writtin in c (not c++). I am using the MinGW-w64 compiler installed via MATLAB add-on. I am working in a 64-bit windows machine with a 64-bit version of MATLAB 2016a.

Thanks, Tim

[Jan]

Please post the code you use to "fill the array coeff". This is the part, which fails, if I understand your question. Then we have to see the corresponding code to get an idea, what's going on.

[Tim Jensen]

The "stripped" version given above also fails. So simply filling the first element by 1.0 causes MATLAB to crash.

Answer 1

[Jan]

What is the type of prhs[8]? Do not rely on "int" to have 32 bits. Better use int32_T oder int64_T specifically.

[Tim Jensen]

I am calling the mex function in MATLAB as:

val = spcomp(uint8(lcomp),lat,lon,h_ell,C,S,GM,a,int32(nmax),int32(mmax));

where the c-file i called "spcomp.c". So prhs[8] and prhs[9] should have 32 bits.