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 in consecutive runs 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,
Panos
void BladeElementGC(const mwIndex InterestBE, const mwIndex BladeElement,
const double *WakeXYZ, const double *RQuarter, const double *BE_RControl, const double *BE_Gamma,
const mwSize NoOfBlades, const mwSize NoOfFilaments, const mwSize NoOfSegments, const mwSize NoOfWakeItems,
const mwSize NoOfElements,
const double conf_CompressibleBiot, const double conf_M, const double conf_CompressibleWakeMachMax,
const double conf_VortexCoreOffset,
const double conf_VortexDelta,
const double conf_kin_visc,
const double conf_VatistasN,
const double conf_BoundVorticity,
double *GCBladeElement)
{
mwIndex Filament;
double Temp1[3], Temp2[3], Temp3[3];
GCBladeElement[0] = 0.;
GCBladeElement[1] = 0.;
GCBladeElement[2] = 0.;
Filament = BladeElement;
TrailingFilamentGC(InterestBE, Filament,
WakeXYZ, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp1);
TrailingFilamentGC(InterestBE, Filament+1,
WakeXYZ, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp2);
if (conf_BoundVorticity>0)
{
BoundElementGC(InterestBE, BladeElement,
WakeXYZ, RQuarter, BE_RControl, BE_Gamma,
NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems,
NoOfElements,
conf_CompressibleBiot, conf_M, conf_CompressibleWakeMachMax,
conf_VortexCoreOffset,
conf_VortexDelta,
conf_kin_visc,
conf_VatistasN,
Temp3);
}
else
{
Temp3[0] = 0.;
Temp3[1] = 0.;
Temp3[2] = 0.;
}
GCBladeElement[0] = Temp2[0] - Temp1[0] + Temp3[0];
GCBladeElement[1] = Temp2[1] - Temp1[1] + Temp3[1];
GCBladeElement[2] = Temp2[2] - Temp1[2] + Temp3[2];
}
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
mwIndex InterestBE;
mwIndex BladeElement;
double *WakeXYZ;
double *RQuarter;
double *BE_RControl;
double *BE_Gamma;
double *GCBladeElement;
const mwSize *dimsWakeXYZ;
const mwSize *dimsBE_RControl;
mwSize NoOfBlades, NoOfFilaments, NoOfSegments, NoOfWakeItems;
mwSize NoOfElements, NoOfDimensions;
double p_CompressibleBiot;
double *p_M;
double *p_CompressibleWakeMachMax;
mwSize p_WakePointsNo;
double 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 = (mwIndex)mxGetScalar(prhs[0]);
/* Get the value of the scalar input Filament */
BladeElement = (mwIndex)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 = (mwSize)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);
}