MEX-fil crasher - C / Matlab problem

Hej Eksperter,

Jeg sidder og er ved at lave et program i MATLAB og jeg har en del af mit program der er seriøst langsomt i forhold til resten. Derfor vil jeg prøve at kode denne i en såkaldt MEX-fil (skrevet i C) for at få det til at køre hurtigere.

Jeg har aldrig før skrevet seriøst i C, men jeg har da fået strikket programmet sammen så det kan compiles. Problemet er bare at hvis jeg herefter kalder MEX-filen et par gange crasher den hele MATLAB. Hvor mange gange svinger lidt men alt fra mellem 5-20 kørsler, og så har vi balladen...

Er der ikke nogen derude der kan se om der er noget hukommelses-vrøvl i min kode eller andet der kan få MATLAB til at crashe.

På forhånd tak,
Killer_R

****************** KODEN KOMMER HER ******************

/*=================================================================
*
* model.c    .MEX file corresponding to model.m
*
* The calling syntax is:
*
*        [S] = model(t, gridp)
*
*=================================================================*/
#include <math.h>
#include "mex.h"

/* Input Arguments */
#define    T_IN    prhs[0]
#define    Y_IN    prhs[1]

/* Output Arguments */
#define    S_OUT    plhs[0]

#if !defined(MAX)
#define    MAX(A, B)    ((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define    MIN(A, B)    ((A) < (B) ? (A) : (B))
#endif

static int len;

static int search(double t[], double value, int low, int high) {
    // Performs a binary search and finds the
    // indicies to that match. Searches indicies
    // between low and high.
    int p;
   
    if (value < t[0]) {
        return -1;
    }
    if (value > t[len]) {
        return len;
    }

    while (low <= high) {
        p = ceil((low+high)/2);
        if (t[p] > value) {
            high = p - 1;
        } else if (t[p] < value) {
            if (value < t[p+1]) {
                return p;
            }
            low = p + 1;
        } else {
            return p;
        }
    }
    return p;
}

static void model(double S[], double gridp[], double t[]) {
    // Extract information:
    double Width = gridp[0];
    double FirstTransit = gridp[1];
    double Period = gridp[2];
    double Depth = gridp[3];
    double tstart = t[0];
    double tend = t[len];
    int j, i, startv, slutv, w, nrtr;
   
    // Udregn tidspunkter for midt-transit:
    double M;
    int nmin = 0;
    int nr = floor((tend-tstart-FirstTransit)/Period) + 1;
    double Mend = tstart+FirstTransit+Period*nr;
    double Msta = tstart+FirstTransit;
    double t1end = Mend - Width/2;
    double t2sta = Msta + Width/2;
    // Dynamiske arrays:
    double *t1, *t2;
    int *k1, *k2;
   
    //mexPrintf("Width: %g \nFirstTransit: %g \nPeriod: %g \n",Width,FirstTransit,Period);
    //mexPrintf("tend: %g \ntstart: %g\n",tend,tstart);
   
    // Tjek af halve transits i starten og slutningen:
    if (Mend+Period > tend && t1end+Period < tend) { // Slutningen
        nr++;
    }
    if (Msta-Period < tstart && t2sta-Period > tstart) { // Starten
        nmin = -1;
    }
    // Giv t1 og t2 deres størrelser:
    nrtr = nr-nmin;
    //mexPrintf("nr: %d\nnmin: %d\nnrtr: %d\n",nr,nmin,nrtr);
    t1 = calloc(nrtr, sizeof(double));
    t2 = calloc(nrtr, sizeof(double));
    k1 = calloc(nrtr, sizeof(int));
    k2 = calloc(nrtr, sizeof(int));
   
    // Find tidspunkterne for transit:
    for (j=nmin; j<=nr; j++) {
        M = tstart+FirstTransit+Period*j;
        t1[j-nmin] = M - Width/2;
        t2[j-nmin] = M + Width/2;
    }
   
    // Den nye metode:
    startv = 0;
    slutv = len;
    for (j=1; j<=nrtr; j++) {
        if (j%2 == 0) {
            w = nrtr+1-j/2;
        } else {
            w = j-floor(j/2);
        }
        w = w-1; // vi er jo 0-baseret!

        k1[w] = search(t, t1[w], startv, slutv) + 1;
        k2[w] = search(t, t2[w], k1[w]+1, slutv);
       
        //mexPrintf("\nw: %d\n k1: %d\n k2: %d\n",w,k1[w],k2[w]);
       
        if (j%2 == 0) {
            slutv = k1[w]-1;
        } else {
            startv = k2[w]+1;
        }
    }
    free(t1);
    free(t2);
   
    // Lav model udfra index:
    for (j=0; j<nrtr; j++) {
        //mexPrintf("j: %d\n",j);
        for (i=k1[j]; i<=k2[j]; i++) {
            //mexPrintf("  i: %d\n",i);
            S[i] = -Depth;
        }
    }
   
    // Clean up:
    free(k1);
    free(k2);
    return;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) {
    double *S, *t, *gridp;
    mwSize m,n;
   
    /* Check for proper number of arguments */
    if (nrhs != 2) {
        mexErrMsgTxt("Two input arguments required.");
    } else if (nlhs > 1) {
        mexErrMsgTxt("Too many output arguments.");
    }
   
    /*  the dimensions of Y.  Y can be 4 X 1 or 1 X 4. */
    m = mxGetM(Y_IN);
    n = mxGetN(Y_IN);
    if (!mxIsDouble(Y_IN) || mxIsComplex(Y_IN) || (MAX(m,n) != 4) || (MIN(m,n) != 1)) {
        mexErrMsgTxt("model requires that gridp be a 4 x 1 vector.");
    }
   
    // Check T:
    m = mxGetM(T_IN);
    n = mxGetN(T_IN);
    if (!mxIsDouble(T_IN) || mxIsComplex(T_IN) || !(MAX(m,n) > 1) || (MIN(m,n) != 1)) {
        mexErrMsgTxt("model requires that t be a n x 1 vector.");
    }
    len = MAX(m,n)-1;
   
    /* Create a matrix for the return argument */
    S_OUT = mxCreateDoubleMatrix(MAX(m,n), 1, mxREAL);
   
    /* Assign pointers to the various parameters */
    S = mxGetPr(S_OUT);
    t = mxGetPr(T_IN);
    gridp = mxGetPr(Y_IN);
       
    /* Do the actual computations in a subroutine */
    model(S, gridp, t);
    return;
}