c - Matlabmexメモリ処理-同じ入力の異なる結果

Question

Matlabコードの計算コストの高い部分をCmexファイルに変換しました。これにより、計算が大幅に高速化されました。問題は、一部の実行ではC mexコードがMatlabコードとまったく同じ結果を（マシンの精度内で）提供するのに対し、他の実行では2つの間に重大な不一致がある（1e-3に達する）ことです。さらに、C mexコードは、まったく同じ入力で異なる結果を連続して実行することにより、それ自体と一致しない場合があります。この問題は、プログラミングによる誤ったメモリ処理が原因である可能性があります。これは、Cmexファイルを作成する最初の試みである可能性が高いです。

計算を実行するために（同じファイル内の）別のC関数を呼び出すmexFunctionをアタッチします。私の問題についての洞察をいただければ幸いです。コードの他の部分が必要な場合はお知らせください。

前もって感謝します、

パノス

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);

 }