optim_private_quad1d.c

ASUFIT－Matlab－全局拟合程序

/*
 * MATLAB Compiler: 2.0.1
 * Date: Tue May 08 21:28:23 2001
 * Arguments: "-B" "sgl" "-m" "-W" "mainhg" "-L" "C" "asufit.m" "absfun.m"
 * "absfunfree.m" "absfunshift.m" "absfunwidth.m" "dispfit.m" "dispfitdisp.m"
 * "dispfitfree.m" "dispfitwidth.m" "seqmodfree.m" "spcfun.m" "spcfun1.m"
 * "atamult.m" "aprecon.m" 
 */
#include "optim_private_quad1d.h"

/*
 * The function "Moptim_private_quad1d" is the implementation version of the
 * "optim/private/quad1d" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\quad1d.m" (lines 1-34). It contains the
 * actual compiled code for that M-function. It is a static function and must
 * only be called from one of the interface functions, appearing below.
 */
/*
 * function[nx,tau] = quad1d(x,ss,delta)
 */
static mxArray * Moptim_private_quad1d(mxArray * * tau,
                                       int nargout_,
                                       mxArray * x,
                                       mxArray * ss,
                                       mxArray * delta) {
    mxArray * nx = mclGetUninitializedArray();
    mxArray * a = mclGetUninitializedArray();
    mxArray * b = mclGetUninitializedArray();
    mxArray * c = mclGetUninitializedArray();
    mxArray * numer = mclGetUninitializedArray();
    mxArray * r1 = mclGetUninitializedArray();
    mxArray * r2 = mclGetUninitializedArray();
    mclValidateInputs("optim/private/quad1d", 3, &x, &ss, &delta);
    /*
     * %QUAD1D	1D quad. zero finder
     * %
     * % [nx,tau] = quad1d(x,ss,delta) tau is min(1,step to zero
     * % of a 1-D quadratic ay^2 + b*y + c , where 
     * % a = x'*x; b = 2*(ss'*x); c = ss'*ss-delta^2). nx is the
     * % new x value, nx = tau*x;
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.2 $  $Date: 1998/03/28 19:59:40 $
     * 
     * a = x'*x;
     */
    mlfAssign(&a, mlfMtimes(mlfCtranspose(x), x));
    /*
     * b = 2*(ss'*x); 
     */
    mlfAssign(&b, mlfMtimes(mlfScalar(2.0), mlfMtimes(mlfCtranspose(ss), x)));
    /*
     * c = ss'*ss-delta^2;
     */
    mlfAssign(
      &c,
      mlfMinus(
        mlfMtimes(mlfCtranspose(ss), ss), mlfMpower(delta, mlfScalar(2.0))));
    /*
     * numer = -(b + sign(b)*sqrt(b^2-4*a*c));
     */
    mlfAssign(
      &numer,
      mlfUminus(
        mlfPlus(
          b,
          mlfMtimes(
            mlfSign(b),
            mlfSqrt(
              mlfMinus(
                mlfMpower(b, mlfScalar(2.0)),
                mlfMtimes(mlfMtimes(mlfScalar(4.0), a), c)))))));
    /*
     * if a == 0 & (numer == 0 | isnan(numer))
     */
    {
        mxArray * a_ = mclInitialize(mlfEq(a, mlfScalar(0.0)));
        if (mlfTobool(a_)) {
            mxArray * b_ = mclInitialize(mlfEq(numer, mlfScalar(0.0)));
            if (mlfTobool(b_)) {
                mlfAssign(&b_, mlfScalar(1));
            } else {
                mlfAssign(&b_, mlfOr(b_, mlfIsnan(numer)));
            }
            {
                int c_0 = mlfTobool(mlfAnd(a_, b_));
                mxDestroyArray(b_);
                if (c_0) {
                    mxDestroyArray(a_);
                    /*
                     * r1 = NaN;
                     */
                    mlfAssign(&r1, mlfNan());
                    goto done_;
                }
            }
        }
        mxDestroyArray(a_);
        /*
         * elseif a == 0 
         */
        if (mlfTobool(mlfEq(a, mlfScalar(0.0)))) {
            /*
             * r1 = inf;
             */
            mlfAssign(&r1, mlfInf());
        /*
         * 
         * else
         */
        } else {
            /*
             * r1 = numer/(2*a);
             */
            mlfAssign(&r1, mlfMrdivide(numer, mlfMtimes(mlfScalar(2.0), a)));
        }
        /*
         * end
         */
        done_:;
    }
    /*
     * r2 = c/(a*r1);
     */
    mlfAssign(&r2, mlfMrdivide(c, mlfMtimes(a, r1)));
    /*
     * tau = max(r1,r2); 
     */
    mlfAssign(tau, mlfMax(NULL, r1, r2, NULL));
    /*
     * tau = min(1,tau);
     */
    mlfAssign(tau, mlfMin(NULL, mlfScalar(1.0), *tau, NULL));
    /*
     * if tau <= 0, 
     */
    if (mlfTobool(mlfLe(*tau, mlfScalar(0.0)))) {
        /*
         * error('square root error in function deljust'); 
         */
        mlfError(mxCreateString("square root error in function deljust"));
    /*
     * end
     */
    }
    /*
     * nx = tau*x;
     */
    mlfAssign(&nx, mlfMtimes(*tau, x));
    mclValidateOutputs("optim/private/quad1d", 2, nargout_, &nx, tau);
    mxDestroyArray(a);
    mxDestroyArray(b);
    mxDestroyArray(c);
    mxDestroyArray(numer);
    mxDestroyArray(r1);
    mxDestroyArray(r2);
    /*
     * 
     * 
     * 
     */
    return nx;
}

/*
 * The function "mlfOptim_private_quad1d" contains the normal interface for the
 * "optim/private/quad1d" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\quad1d.m" (lines 1-34). This function
 * processes any input arguments and passes them to the implementation version
 * of the function, appearing above.
 */
mxArray * mlfOptim_private_quad1d(mxArray * * tau,
                                  mxArray * x,
                                  mxArray * ss,
                                  mxArray * delta) {
    int nargout = 1;
    mxArray * nx = mclGetUninitializedArray();
    mxArray * tau__ = mclGetUninitializedArray();
    mlfEnterNewContext(1, 3, tau, x, ss, delta);
    if (tau != NULL) {
        ++nargout;
    }
    nx = Moptim_private_quad1d(&tau__, nargout, x, ss, delta);
    mlfRestorePreviousContext(1, 3, tau, x, ss, delta);
    if (tau != NULL) {
        mclCopyOutputArg(tau, tau__);
    } else {
        mxDestroyArray(tau__);
    }
    return mlfReturnValue(nx);
}

/*
 * The function "mlxOptim_private_quad1d" contains the feval interface for the
 * "optim/private/quad1d" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\quad1d.m" (lines 1-34). The feval
 * function calls the implementation version of optim/private/quad1d through
 * this function. This function processes any input arguments and passes them
 * to the implementation version of the function, appearing above.
 */
void mlxOptim_private_quad1d(int nlhs,
                             mxArray * plhs[],
                             int nrhs,
                             mxArray * prhs[]) {
    mxArray * mprhs[3];
    mxArray * mplhs[2];
    int i;
    if (nlhs > 2) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: optim/private/quad1d Line: 1 Col"
            "umn: 0 The function \"optim/private/quad1d\" was calle"
            "d with more than the declared number of outputs (2)"));
    }
    if (nrhs > 3) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: optim/private/quad1d Line: 1 Col"
            "umn: 0 The function \"optim/private/quad1d\" was calle"
            "d with more than the declared number of inputs (3)"));
    }
    for (i = 0; i < 2; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 3 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 3; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(0, 3, mprhs[0], mprhs[1], mprhs[2]);
    mplhs[0]
      = Moptim_private_quad1d(&mplhs[1], nlhs, mprhs[0], mprhs[1], mprhs[2]);
    mlfRestorePreviousContext(0, 3, mprhs[0], mprhs[1], mprhs[2]);
    plhs[0] = mplhs[0];
    for (i = 1; i < 2 && i < nlhs; ++i) {
        plhs[i] = mplhs[i];
    }
    for (; i < 2; ++i) {
        mxDestroyArray(mplhs[i]);
    }
}