optim_private_trust.c

ASUFIT－Matlab－全局拟合程序

    mclValidateInputs("trust/seceqn", 4, &lambda, &eigval, &alpha, &delta);
    /*
     * %SEC	Secular equation
     * %
     * % value = SEC(lambda,eigval,alpha,delta) returns the value
     * % of the secular equation at a set of m points lambda
     * %
     * %
     * 
     * 
     * %
     * m = length(lambda); n = length(eigval);
     */
    mlfAssign(&m, mlfLength(lambda));
    mlfAssign(&n, mlfLength(eigval));
    /*
     * unn = ones(n,1); unm = ones(m,1);
     */
    mlfAssign(&unn, mlfOnes(n, mlfScalar(1.0), NULL));
    mlfAssign(&unm, mlfOnes(m, mlfScalar(1.0), NULL));
    /*
     * M = eigval*unm' + unn*lambda'; MC = M;
     */
    mlfAssign(
      &M,
      mlfPlus(
        mlfMtimes(eigval, mlfCtranspose(unm)),
        mlfMtimes(unn, mlfCtranspose(lambda))));
    mlfAssign(&MC, M);
    /*
     * MM = alpha*unm';
     */
    mlfAssign(&MM, mlfMtimes(alpha, mlfCtranspose(unm)));
    /*
     * M(M~=0) = MM(M~=0) ./ M(M~=0);
     */
    mlfIndexAssign(
      &M,
      "(?)",
      mlfNe(M, mlfScalar(0.0)),
      mlfRdivide(
        mlfIndexRef(MM, "(?)", mlfNe(M, mlfScalar(0.0))),
        mlfIndexRef(M, "(?)", mlfNe(M, mlfScalar(0.0)))));
    /*
     * M(MC==0) = Inf*ones(size(MC(MC==0))); 
     */
    mlfIndexAssign(
      &M,
      "(?)",
      mlfEq(MC, mlfScalar(0.0)),
      mlfMtimes(
        mlfInf(),
        mlfOnes(
          mlfSize(
            mclValueVarargout(),
            mlfIndexRef(MC, "(?)", mlfEq(MC, mlfScalar(0.0))),
            NULL),
          NULL)));
    /*
     * M = M.*M;
     */
    mlfAssign(&M, mlfTimes(M, M));
    /*
     * value = sqrt(unm ./ (M'*unn));
     */
    mlfAssign(
      &value, mlfSqrt(mlfRdivide(unm, mlfMtimes(mlfCtranspose(M), unn))));
    /*
     * value(value==NaN) = zeros(length(value(value==NaN)),1);
     */
    mlfIndexAssign(
      &value,
      "(?)",
      mlfEq(value, mlfNan()),
      mlfZeros(
        mlfLength(mlfIndexRef(value, "(?)", mlfEq(value, mlfNan()))),
        mlfScalar(1.0),
        NULL));
    /*
     * value = (1/delta)*unm - value;
     */
    mlfAssign(
      &value,
      mlfMinus(mlfMtimes(mlfMrdivide(mlfScalar(1.0), delta), unm), value));
    mclValidateOutputs("trust/seceqn", 1, nargout_, &value);
    mxDestroyArray(M);
    mxDestroyArray(MC);
    mxDestroyArray(MM);
    mxDestroyArray(m);
    mxDestroyArray(n);
    mxDestroyArray(unm);
    mxDestroyArray(unn);
    /*
     * 
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * 
     */
    return value;
}

/*
 * The function "mlfTrust_seceqn" contains the normal interface for the
 * "trust/seceqn" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 114-138). This function
 * processes any input arguments and passes them to the implementation version
 * of the function, appearing above.
 */
static mxArray * mlfTrust_seceqn(mxArray * lambda,
                                 mxArray * eigval,
                                 mxArray * alpha,
                                 mxArray * delta) {
    int nargout = 1;
    mxArray * value = mclGetUninitializedArray();
    mlfEnterNewContext(0, 4, lambda, eigval, alpha, delta);
    value = Mtrust_seceqn(nargout, lambda, eigval, alpha, delta);
    mlfRestorePreviousContext(0, 4, lambda, eigval, alpha, delta);
    return mlfReturnValue(value);
}

/*
 * The function "mlxTrust_seceqn" contains the feval interface for the
 * "trust/seceqn" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 114-138). The feval
 * function calls the implementation version of trust/seceqn through this
 * function. This function processes any input arguments and passes them to the
 * implementation version of the function, appearing above.
 */
static void mlxTrust_seceqn(int nlhs,
                            mxArray * plhs[],
                            int nrhs,
                            mxArray * prhs[]) {
    mxArray * mprhs[4];
    mxArray * mplhs[1];
    int i;
    if (nlhs > 1) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: trust/seceqn Line: 114 Colum"
            "n: 0 The function \"trust/seceqn\" was called with"
            " more than the declared number of outputs (1)"));
    }
    if (nrhs > 4) {
        mlfError(
          mxCreateString(
            "Run-time Error: File: trust/seceqn Line: 114 Colu"
            "mn: 0 The function \"trust/seceqn\" was called wi"
            "th more than the declared number of inputs (4)"));
    }
    for (i = 0; i < 1; ++i) {
        mplhs[i] = NULL;
    }
    for (i = 0; i < 4 && i < nrhs; ++i) {
        mprhs[i] = prhs[i];
    }
    for (; i < 4; ++i) {
        mprhs[i] = NULL;
    }
    mlfEnterNewContext(0, 4, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
    mplhs[0] = Mtrust_seceqn(nlhs, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
    mlfRestorePreviousContext(0, 4, mprhs[0], mprhs[1], mprhs[2], mprhs[3]);
    plhs[0] = mplhs[0];
}

/*
 * The function "Mtrust_rfzero" is the implementation version of the
 * "trust/rfzero" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 138-290). 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 [b,c,itfun] = rfzero(FunFcn,x,itbnd,eigval,alpha,delta,tol,trace)
 */
static mxArray * Mtrust_rfzero(mxArray * * c,
                               mxArray * * itfun,
                               int nargout_,
                               mxArray * FunFcn,
                               mxArray * x,
                               mxArray * itbnd,
                               mxArray * eigval,
                               mxArray * alpha,
                               mxArray * delta,
                               mxArray * tol,
                               mxArray * trace) {
    mxArray * b = mclGetUninitializedArray();
    mxArray * a = mclGetUninitializedArray();
    mxArray * d = mclGetUninitializedArray();
    mxArray * dx = mclGetUninitializedArray();
    mxArray * e = mclGetUninitializedArray();
    mxArray * fa = mclGetUninitializedArray();
    mxArray * fb = mclGetUninitializedArray();
    mxArray * fc = mclGetUninitializedArray();
    mxArray * init = mclGetUninitializedArray();
    mxArray * m = mclGetUninitializedArray();
    mxArray * nargin_ = mclGetUninitializedArray();
    mxArray * p = mclGetUninitializedArray();
    mxArray * q = mclGetUninitializedArray();
    mxArray * r = mclGetUninitializedArray();
    mxArray * s = mclGetUninitializedArray();
    mxArray * sign = mclGetUninitializedArray();
    mxArray * step = mclGetUninitializedArray();
    mxArray * toler = mclGetUninitializedArray();
    mlfAssign(
      &nargin_,
      mlfNargin(0, FunFcn, x, itbnd, eigval, alpha, delta, tol, trace, NULL));
    mclValidateInputs(
      "trust/rfzero",
      8,
      &FunFcn,
      &x,
      &itbnd,
      &eigval,
      &alpha,
      &delta,
      &tol,
      &trace);
    mclCopyArray(&tol);
    mclCopyArray(&trace);
    /*
     * %RFZERO Find zero to the right
     * %
     * %	[b,c,itfun] = rfzero(FunFcn,x,itbnd,eigval,alpha,delta,tol,trace)
     * %	Zero of a function of one variable to the RIGHT of the
     * %       starting point x. A small modification of the M-file fzero,
     * %       described below, to ensure a zero to the Right of x is
     * %       searched for.
     * %
     * %	RFZERO is a slightly modified version of function FZERO
     * 
     * 
     * 
     * 
     * %	FZERO(F,X) finds a zero of f(x).  F is a string containing the
     * %	name of a real-valued function of a single real variable.   X is
     * %	a starting guess.  The value returned is near a point where F
     * %	changes sign.  For example, FZERO('sin',3) is pi.  Note the
     * %	quotes around sin.  Ordinarily, functions are defined in M-files.
     * %
     * %	An optional third argument sets the relative tolerance for the
     * %	convergence test.   The presence of an nonzero optional fourth
     * %	argument triggers a printing trace of the steps.
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %	C.B. Moler 1-19-86
     * %	Revised CBM 3-25-87, LS 12-01-88.
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * %  This algorithm was originated by T. Dekker.  An Algol 60 version,
     * %  with some improvements, is given by Richard Brent in "Algorithms for
     * %  Minimization Without Derivatives", Prentice-Hall, 1973.  A Fortran
     * %  version is in Forsythe, Malcolm and Moler, "Computer Methods
     * %  for Mathematical Computations", Prentice-Hall, 1976.
     * %
     * % Initialization
     * if nargin < 7, trace = 0; tol = eps; end
     */
    if (mlfTobool(mlfLt(nargin_, mlfScalar(7.0)))) {
        mlfAssign(&trace, mlfScalar(0.0));
        mlfAssign(&tol, mlfEps());
    }
    /*
     * if nargin == 7, trace = 0; end
     */
    if (mlfTobool(mlfEq(nargin_, mlfScalar(7.0)))) {
        mlfAssign(&trace, mlfScalar(0.0));
    }
    /*
     * if trace, clc, end
     */
    if (mlfTobool(trace)) {
        mlfClc();
    }
    /*
     * itfun = 0;
     */
    mlfAssign(itfun, mlfScalar(0.0));
    /*
     * %
     * %if x ~= 0, dx = x/20;
     * %if x ~= 0, dx = abs(x)/20;
     * if x~= 0, dx = abs(x)/2;
     */
    if (mlfTobool(mlfNe(x, mlfScalar(0.0)))) {
        mlfAssign(&dx, mlfMrdivide(mlfAbs(x), mlfScalar(2.0)));
    /*
     * %
     * %else, dx = 1/20;
     * else, dx = 1/2;
     */
    } else {
        mlfAssign(&dx, mlfMrdivide(mlfScalar(1.0), mlfScalar(2.0)));
    /*
     * end
     */
    }
    /*
     * %
     * %a = x - dx;  fa = feval(FunFcn,a,eigval,alpha,delta);
     * a = x; c = a;  fa = feval(FunFcn,a,eigval,alpha,delta);
     */
    mlfAssign(&a, x);
    mlfAssign(c, a);
    mlfAssign(
      &fa,
      mlfFeval(
        mclValueVarargout(),
        mclFevalLookup(FunFcn, 2, local_function_table_),
        a,
        eigval,
        alpha,
        delta,
        NULL));
    /*
     * itfun = itfun+1;
     */
    mlfAssign(itfun, mlfPlus(*itfun, mlfScalar(1.0)));
    /*
     * %
     * if trace
     */
    if (mlfTobool(trace)) {
        /*
         * %home
         * init = [a fa]
         */
        mlfAssign(&init, mlfHorzcat(a, fa, NULL));
        mclPrintArray(init, "init");
    /*
     * end
     */
    }
    /*
     * b = x + dx;
     */
    mlfAssign(&b, mlfPlus(x, dx));
    /*
     * b = x + 1;
     */
    mlfAssign(&b, mlfPlus(x, mlfScalar(1.0)));
    /*
     * fb = feval(FunFcn,b,eigval,alpha,delta);
     */
    mlfAssign(
      &fb,
      mlfFeval(
        mclValueVarargout(),
        mclFevalLookup(FunFcn, 2, local_function_table_),
        b,
        eigval,
        alpha,
        delta,
        NULL));
    /*
     * itfun = itfun+1;
     */
    mlfAssign(itfun, mlfPlus(*itfun, mlfScalar(1.0)));
    /*
     * if trace
     */
    if (mlfTobool(trace)) {
        /*
         * %home
         * init = [b fb]
         */
        mlfAssign(&init, mlfHorzcat(b, fb, NULL));
        mclPrintArray(init, "init");
    /*
     * end
     */
    }
    /*
     * 
     * % Find change of sign.
     * 
     * while (fa > 0) == (fb > 0)
     */
    while (mlfTobool(
             mlfEq(mlfGt(fa, mlfScalar(0.0)), mlfGt(fb, mlfScalar(0.0))))) {
        /*
         * dx = 2*dx;
         */
        mlfAssign(&dx, mlfMtimes(mlfScalar(2.0), dx));
        /*
         * %
         * %  a = x - dx;  fa = feval(FunFcn,a);
         * %  if trace, home, sign = [a fa], end
         * %
         * if (fa > 0) ~= (fb > 0), break, end
         */
        if (mlfTobool(
              mlfNe(mlfGt(fa, mlfScalar(0.0)), mlfGt(fb, mlfScalar(0.0))))) {
            break;
        }
        /*
         * b = x + dx;  fb = feval(FunFcn,b,eigval,alpha,delta);
         */
        mlfAssign(&b, mlfPlus(x, dx));
        mlfAssign(
          &fb,
          mlfFeval(
            mclValueVarargout(),
            mclFevalLookup(FunFcn, 2, local_function_table_),
            b,
            eigval,
            alpha,
            delta,
            NULL));
        /*
         * itfun = itfun+1;
         */
        mlfAssign(itfun, mlfPlus(*itfun, mlfScalar(1.0)));
        /*
         * if trace
         */
        if (mlfTobool(trace)) {
            /*
             * %home
             * sign = [b fb]
             */
            mlfAssign(&sign, mlfHorzcat(b, fb, NULL));
            mclPrintArray(sign, "sign");
        /*
         * end
         */
        }
        /*
         * if itfun > itbnd, break; end
         */
        if (mlfTobool(mlfGt(*itfun, itbnd))) {
            break;
        }
    /*
     * end
     */
    }
    /*
     *