Hello all,

I ve converted a computationally expensive part of my Matlab code into a C mex file, which greatly speeded up the calculation. The problem is that in some runs the C mex code gives exactly the same result as the Matlab code (within the machine accuracy), whereas in some others there is a significant discrepancy between the two (reaching even 1e-3). Furthermore, there are times that the C mex code disagrees with itself by giving different results with exactly the same inputs. I guess that the problem is probably caused by some incorrect memory handling due to my programming. This is highly likely as thats my first attempt to create a C mex file.

I attach my mexFunction which calls another C function (in the same file) to perform the calculations. I would be grateful if you could give me some insight on my problem. Please let me know if you need any other part of the code.

Thanks in advance
Panagiotis

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    
    int InterestBE;
    int BladeElement;
    double *WakeXYZ;
    double *RQuarter;
    double *BE_RControl;
    double *BE_Gamma;
    const mxArray *Conf;
    double *GCBladeElement;
   
    const int *dimsWakeXYZ;
    const int *dimsBE_RControl;
    
    int NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems;
    int NoOfElements, NoOfDimensions;
    
    int p_CompressibleBiot;
    double *p_M;
    double *p_CompressibleWakeMachMax;
    int p_WakePointsNo;
    int p_BoundVorticity;
    
    double *p_VortexCoreOffset;
    double *p_VortexDelta;
    double *p_kin_visc;
    double *p_VatistasN;
    
    
    
    //----------------------------------------------------------
    // First retrieve inputs
    //----------------------------------------------------------
    
    /* Get the value of the scalar input InterestBE  */
    InterestBE = mxGetScalar(prhs[0]);
    
    /* Get the value of the scalar input Filament  */
    BladeElement = mxGetScalar(prhs[1]);
    
    /* create a pointer to the real data in the WakeXYZ matrix  */
    WakeXYZ = mxGetPr(prhs[2]);
    
    /* create a pointer to the real data in the RQuarter matrix  */
    RQuarter = mxGetPr(prhs[3]);
    
    /* create a pointer to the real data in the BE_RControl matrix  */
    BE_RControl = mxGetPr(prhs[4]);
    
    /* create a pointer to the real data in the BE_RControl matrix  */
    BE_Gamma = mxGetPr(prhs[5]);
    
    
    /* create pointers to the fields of the conf structure  */
    
    p_CompressibleBiot = mxGetScalar(mxGetField(prhs[6], 0, "CompressibleBiot"));
    p_M = mxGetPr(mxGetField(prhs[6], 0, "M"));
    p_CompressibleWakeMachMax = mxGetPr(mxGetField(prhs[6], 0, "CompressibleWakeMachMax"));
    p_VortexCoreOffset = mxGetPr(mxGetField(prhs[6], 0, "VortexCoreOffset"));
    p_VortexDelta = mxGetPr(mxGetField(prhs[6], 0, "VortexDelta"));
    p_kin_visc = mxGetPr(mxGetField(prhs[6], 0, "kin_visc"));
    p_VatistasN = mxGetPr(mxGetField(prhs[6], 0, "VatistasN"));
    p_WakePointsNo = mxGetScalar(mxGetField(prhs[6],0,"WakePoints"));
    p_BoundVorticity = mxGetScalar(mxGetField(prhs[6],0,"BoundVorticity"));
    
    
    // Get the dimensions of the WakeXYZ array
    dimsWakeXYZ = mxGetDimensions(prhs[2]);

    // Get the dimensions of the BE_RControl array
    dimsBE_RControl = mxGetDimensions(prhs[4]);
    
    NoOfBlades = dimsWakeXYZ[0];
    NoOfFilaments = dimsWakeXYZ[1];
    NoOfSegments = dimsWakeXYZ[2];
    NoOfWakeItems = dimsWakeXYZ[3];
    NoOfElements = dimsBE_RControl[1];
    
    //----------------------------------------------------------
    // Call the calculation subroutine
    //----------------------------------------------------------
    
    /* create the output matrix */
    plhs[0] = mxCreateDoubleMatrix(1,3,mxREAL);
    GCBladeElement = mxGetPr(plhs[0]);
    
    BladeElementGC(InterestBE, BladeElement,
                WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
                NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
                NoOfElements,
                p_CompressibleBiot, *p_M, *p_CompressibleWakeMachMax,
                *p_VortexCoreOffset,
                *p_VortexDelta,
                *p_kin_visc,
                *p_VatistasN,
                p_BoundVorticity,
                GCBladeElement);
       
}

"Panagiotis" wrote in message <jrhk43$hf2$1@newscl01ah.mathworks.com>...
>     
>     BladeElementGC(InterestBE, BladeElement,
>                 WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
>                 NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
>                 NoOfElements,
>                 p_CompressibleBiot, *p_M, *p_CompressibleWakeMachMax,
>                 *p_VortexCoreOffset,
>                 *p_VortexDelta,
>                 *p_kin_visc,
>                 *p_VatistasN,
>                 p_BoundVorticity,
>                 GCBladeElement);

What does the BladeElementGC interface look like? Normally I would not expect you to be dereferenceing the pointers in the argument list. I.e., why are you passing *p_kin_visc instead of just p_kin_visc? If you are accessing array elements inside the BladeElementGC routine then dereferencing the pointers in the argument list will not work. But again, I would need to see the actual interface to be sure what is going wrong.

James Tursa

Dear James, I have sent my reply on Matlab Answers.

thanks
Panos