optim_private_trust.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_trust.h"

static mxArray * Moptim_private_trust(mxArray * * val,
                                      mxArray * * posdef,
                                      mxArray * * count,
                                      mxArray * * lambda,
                                      int nargout_,
                                      mxArray * g,
                                      mxArray * H,
                                      mxArray * delta);
static mxArray * Mtrust_seceqn(int nargout_,
                               mxArray * lambda,
                               mxArray * eigval,
                               mxArray * alpha,
                               mxArray * delta);
static mxArray * mlfTrust_seceqn(mxArray * lambda,
                                 mxArray * eigval,
                                 mxArray * alpha,
                                 mxArray * delta);
static void mlxTrust_seceqn(int nlhs,
                            mxArray * plhs[],
                            int nrhs,
                            mxArray * prhs[]);
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);
static mxArray * mlfTrust_rfzero(mxArray * * c,
                                 mxArray * * itfun,
                                 mxArray * FunFcn,
                                 mxArray * x,
                                 mxArray * itbnd,
                                 mxArray * eigval,
                                 mxArray * alpha,
                                 mxArray * delta,
                                 mxArray * tol,
                                 mxArray * trace);
static void mlxTrust_rfzero(int nlhs,
                            mxArray * plhs[],
                            int nrhs,
                            mxArray * prhs[]);

static mlfFunctionTableEntry local_function_table_[2]
  = { { "seceqn", mlxTrust_seceqn, 4, 1 },
      { "rfzero", mlxTrust_rfzero, 8, 3 } };

/*
 * The function "Moptim_private_trust" is the implementation version of the
 * "optim/private/trust" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\trust.m" (lines 1-114). 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 [s,val,posdef,count,lambda] = trust(g,H,delta)
 */
static mxArray * Moptim_private_trust(mxArray * * val,
                                      mxArray * * posdef,
                                      mxArray * * count,
                                      mxArray * * lambda,
                                      int nargout_,
                                      mxArray * g,
                                      mxArray * H,
                                      mxArray * delta) {
    mxArray * s = mclGetUninitializedArray();
    mxArray * D = mclGetUninitializedArray();
    mxArray * V = mclGetUninitializedArray();
    mxArray * alpha = mclGetUninitializedArray();
    mxArray * arg = mclGetUninitializedArray();
    mxArray * arg1 = mclGetUninitializedArray();
    mxArray * arg2 = mclGetUninitializedArray();
    mxArray * b = mclGetUninitializedArray();
    mxArray * beta = mclGetUninitializedArray();
    mxArray * c = mclGetUninitializedArray();
    mxArray * coeff = mclGetUninitializedArray();
    mxArray * eigval = mclGetUninitializedArray();
    mxArray * itbnd = mclGetUninitializedArray();
    mxArray * jmin = mclGetUninitializedArray();
    mxArray * key = mclGetUninitializedArray();
    mxArray * lam = mclGetUninitializedArray();
    mxArray * laminit = mclGetUninitializedArray();
    mxArray * mineig = mclGetUninitializedArray();
    mxArray * n = mclGetUninitializedArray();
    mxArray * nrms = mclGetUninitializedArray();
    mxArray * sig = mclGetUninitializedArray();
    mxArray * tol = mclGetUninitializedArray();
    mxArray * tol2 = mclGetUninitializedArray();
    mxArray * vval = mclGetUninitializedArray();
    mxArray * w = mclGetUninitializedArray();
    mclValidateInputs("optim/private/trust", 3, &g, &H, &delta);
    mclCopyArray(&H);
    /*
     * %TRUST	Exact soln of trust region problem
     * %
     * % [s,val,posdef,count,lambda] = TRUST(g,H,delta) Solves the trust region
     * % problem: min{g^Ts + 1/2 s^THs: ||s|| <= delta}. The full
     * % eigen-decomposition is used; based on the secular equation,
     * % 1/delta - 1/(||s||) = 0. The solution is s, the value
     * % of the quadratic at the solution is val; posdef = 1 if H
     * % is pos. definite; otherwise posdef = 0. The number of
     * % evaluations of the secular equation is count, lambda
     * % is the value of the corresponding Lagrange multiplier.
     * %
     * %
     * % TRUST is meant to be applied to very small dimensional problems.
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.3 $  $Date: 1998/10/09 14:02:49 $
     * 
     * % INITIALIZATION
     * tol = 10^(-12); 
     */
    mlfAssign(&tol, mlfMpower(mlfScalar(10.0), mlfScalar(-12.0)));
    /*
     * tol2 = 10^(-8); 
     */
    mlfAssign(&tol2, mlfMpower(mlfScalar(10.0), mlfScalar(-8.0)));
    /*
     * key = 0; 
     */
    mlfAssign(&key, mlfScalar(0.0));
    /*
     * itbnd = 50;
     */
    mlfAssign(&itbnd, mlfScalar(50.0));
    /*
     * lambda = 0;
     */
    mlfAssign(lambda, mlfScalar(0.0));
    /*
     * n = length(g); 
     */
    mlfAssign(&n, mlfLength(g));
    /*
     * coeff(1:n,1) = zeros(n,1);
     */
    mlfIndexAssign(
      &coeff,
      "(?,?)",
      mlfColon(mlfScalar(1.0), n, NULL),
      mlfScalar(1.0),
      mlfZeros(n, mlfScalar(1.0), NULL));
    /*
     * H = full(H);
     */
    mlfAssign(&H, mlfFull(H));
    /*
     * [V,D] = eig(H); 
     */
    mlfAssign(&V, mlfEig(&D, H, NULL));
    /*
     * count = 0;
     */
    mlfAssign(count, mlfScalar(0.0));
    /*
     * eigval = diag(D); 
     */
    mlfAssign(&eigval, mlfDiag(D, NULL));
    /*
     * [mineig,jmin] = min(eigval);
     */
    mlfAssign(&mineig, mlfMin(&jmin, eigval, NULL, NULL));
    /*
     * alpha = -V'*g; 
     */
    mlfAssign(&alpha, mlfMtimes(mlfUminus(mlfCtranspose(V)), g));
    /*
     * sig = sign(alpha(jmin)) + (alpha(jmin)==0);
     */
    mlfAssign(
      &sig,
      mlfPlus(
        mlfSign(mlfIndexRef(alpha, "(?)", jmin)),
        mlfEq(mlfIndexRef(alpha, "(?)", jmin), mlfScalar(0.0))));
    /*
     * 
     * % POSITIVE DEFINITE CASE
     * if mineig > 0
     */
    if (mlfTobool(mlfGt(mineig, mlfScalar(0.0)))) {
        /*
         * coeff = alpha ./ eigval; 
         */
        mlfAssign(&coeff, mlfRdivide(alpha, eigval));
        /*
         * lambda = 0;
         */
        mlfAssign(lambda, mlfScalar(0.0));
        /*
         * s = V*coeff; 
         */
        mlfAssign(&s, mlfMtimes(V, coeff));
        /*
         * posdef = 1; 
         */
        mlfAssign(posdef, mlfScalar(1.0));
        /*
         * nrms = norm(s);
         */
        mlfAssign(&nrms, mlfNorm(s, NULL));
        /*
         * if nrms <= 1.2*delta
         */
        if (mlfTobool(mlfLe(nrms, mlfMtimes(mlfScalar(1.2), delta)))) {
            /*
             * key = 1;
             */
            mlfAssign(&key, mlfScalar(1.0));
        /*
         * else 
         */
        } else {
            /*
             * laminit = 0;
             */
            mlfAssign(&laminit, mlfScalar(0.0));
        /*
         * end
         */
        }
    /*
     * else
     */
    } else {
        /*
         * laminit = -mineig; 
         */
        mlfAssign(&laminit, mlfUminus(mineig));
        /*
         * posdef = 0;
         */
        mlfAssign(posdef, mlfScalar(0.0));
    /*
     * end
     */
    }
    /*
     * 
     * % INDEFINITE CASE
     * if key == 0
     */
    if (mlfTobool(mlfEq(key, mlfScalar(0.0)))) {
        /*
         * if seceqn(laminit,eigval,alpha,delta) > 0
         */
        if (mlfTobool(
              mlfGt(
                mlfTrust_seceqn(laminit, eigval, alpha, delta),
                mlfScalar(0.0)))) {
            /*
             * [b,c,count] = rfzero('seceqn',laminit,itbnd,eigval,alpha,delta,tol);
             */
            mlfAssign(
              &b,
              mlfTrust_rfzero(
                &c,
                count,
                mxCreateString("seceqn"),
                laminit,
                itbnd,
                eigval,
                alpha,
                delta,
                tol,
                NULL));
            /*
             * vval = abs(seceqn(b,eigval,alpha,delta));
             */
            mlfAssign(&vval, mlfAbs(mlfTrust_seceqn(b, eigval, alpha, delta)));
            /*
             * if abs(seceqn(b,eigval,alpha,delta)) <= tol2
             */
            if (mlfTobool(
                  mlfLe(
                    mlfAbs(mlfTrust_seceqn(b, eigval, alpha, delta)), tol2))) {
                /*
                 * lambda = b; 
                 */
                mlfAssign(lambda, b);
                /*
                 * key = 2;
                 */
                mlfAssign(&key, mlfScalar(2.0));
                /*
                 * lam = lambda*ones(n,1);
                 */
                mlfAssign(
                  &lam, mlfMtimes(*lambda, mlfOnes(n, mlfScalar(1.0), NULL)));
                /*
                 * w = eigval + lam;
                 */
                mlfAssign(&w, mlfPlus(eigval, lam));
                /*
                 * arg1 = (w==0) & (alpha == 0); 
                 */
                mlfAssign(
                  &arg1,
                  mlfAnd(
                    mlfEq(w, mlfScalar(0.0)), mlfEq(alpha, mlfScalar(0.0))));
                /*
                 * arg2 = (w==0) & (alpha ~= 0);
                 */
                mlfAssign(
                  &arg2,
                  mlfAnd(
                    mlfEq(w, mlfScalar(0.0)), mlfNe(alpha, mlfScalar(0.0))));
                /*
                 * coeff(w ~=0) = alpha(w ~=0) ./ w(w~=0);
                 */
                mlfIndexAssign(
                  &coeff,
                  "(?)",
                  mlfNe(w, mlfScalar(0.0)),
                  mlfRdivide(
                    mlfIndexRef(alpha, "(?)", mlfNe(w, mlfScalar(0.0))),
                    mlfIndexRef(w, "(?)", mlfNe(w, mlfScalar(0.0)))));
                /*
                 * coeff(arg1) = 0;
                 */
                mlfIndexAssign(&coeff, "(?)", arg1, mlfScalar(0.0));
                /*
                 * coeff(arg2) = Inf;
                 */
                mlfIndexAssign(&coeff, "(?)", arg2, mlfInf());
                /*
                 * coeff(isnan(coeff))=0;
                 */
                mlfIndexAssign(&coeff, "(?)", mlfIsnan(coeff), mlfScalar(0.0));
                /*
                 * s = V*coeff; 
                 */
                mlfAssign(&s, mlfMtimes(V, coeff));
                /*
                 * nrms = norm(s);
                 */
                mlfAssign(&nrms, mlfNorm(s, NULL));
                /*
                 * if (nrms > 1.2*delta) | (nrms < .8*delta)
                 */
                {
                    mxArray * a_ = mclInitialize(
                                     mlfGt(
                                       nrms, mlfMtimes(mlfScalar(1.2), delta)));
                    if (mlfTobool(a_)
                        || mlfTobool(
                             mlfOr(
                               a_,
                               mlfLt(
                                 nrms, mlfMtimes(mlfScalar(.8), delta))))) {
                        mxDestroyArray(a_);
                        /*
                         * key = 5; 
                         */
                        mlfAssign(&key, mlfScalar(5.0));
                        /*
                         * lambda = -mineig;
                         */
                        mlfAssign(lambda, mlfUminus(mineig));
                    } else {
                        mxDestroyArray(a_);
                    }
                /*
                 * end
                 */
                }
            /*
             * else
             */
            } else {
                /*
                 * lambda = -mineig; 
                 */
                mlfAssign(lambda, mlfUminus(mineig));
                /*
                 * key = 3;
                 */
                mlfAssign(&key, mlfScalar(3.0));
            /*
             * end
             */
            }
        /*
         * else
         */
        } else {
            /*
             * lambda = -mineig; 
             */
            mlfAssign(lambda, mlfUminus(mineig));
            /*
             * key = 4;
             */
            mlfAssign(&key, mlfScalar(4.0));
        /*
         * end
         */
        }
        /*
         * lam = lambda*ones(n,1);
         */
        mlfAssign(&lam, mlfMtimes(*lambda, mlfOnes(n, mlfScalar(1.0), NULL)));
        /*
         * if (key > 2) 
         */
        if (mlfTobool(mlfGt(key, mlfScalar(2.0)))) {
            /*
             * arg = abs(eigval + lam) < 10 * eps * max(abs(eigval),ones(n,1));
             */
            mlfAssign(
              &arg,
              mlfLt(
                mlfAbs(mlfPlus(eigval, lam)),
                mlfMtimes(
                  mlfMtimes(mlfScalar(10.0), mlfEps()),
                  mlfMax(
                    NULL,
                    mlfAbs(eigval),
                    mlfOnes(n, mlfScalar(1.0), NULL),
                    NULL))));
            /*
             * alpha(arg) = 0;