optim_private_snls.c

ASUFIT－Matlab－全局拟合程序

/*
 * MATLAB Compiler: 2.0.1
 * Date: Tue May 08 21:28:22 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_snls.h"
#include "color.h"
#include "optim_private_definev.h"
#include "optim_private_display1.h"
#include "optim_private_perturb.h"
#include "sfdnls.h"
#include "optim_private_startx.h"
#include "optim_private_trdog.h"
#include "optim_private_xplot.h"

/*
 * The function "Moptim_private_snls" is the implementation version of the
 * "optim/private/snls" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\snls.m" (lines 1-424). 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[xcurr,fvec,LAMBDA,JACOB,EXITFLAG,OUTPUT]=snls(funfcn,xstart,l,u,verb,options,fval,JACval,YDATA,caller,Jstr,computeLambda,varargin)
 */
static mxArray * Moptim_private_snls(mxArray * * fvec,
                                     mxArray * * LAMBDA,
                                     mxArray * * JACOB,
                                     mxArray * * EXITFLAG,
                                     mxArray * * OUTPUT,
                                     int nargout_,
                                     mxArray * funfcn,
                                     mxArray * xstart,
                                     mxArray * l,
                                     mxArray * u,
                                     mxArray * verb,
                                     mxArray * options,
                                     mxArray * fval,
                                     mxArray * JACval,
                                     mxArray * YDATA,
                                     mxArray * caller,
                                     mxArray * Jstr,
                                     mxArray * computeLambda,
                                     mxArray * varargin) {
    mxArray * xcurr = mclGetUninitializedArray();
    mxArray * A = mclGetUninitializedArray();
    mxArray * D = mclGetUninitializedArray();
    mxArray * DS = mclGetUninitializedArray();
    mxArray * Z = mclGetUninitializedArray();
    mxArray * active_tol = mclGetUninitializedArray();
    mxArray * ans = mclInitializeAns();
    mxArray * arg = mclGetUninitializedArray();
    mxArray * arg2 = mclGetUninitializedArray();
    mxArray * argl = mclGetUninitializedArray();
    mxArray * argu = mclGetUninitializedArray();
    mxArray * aug = mclGetUninitializedArray();
    mxArray * bndfeas = mclGetUninitializedArray();
    mxArray * currOutput = mclGetUninitializedArray();
    mxArray * dd = mclGetUninitializedArray();
    mxArray * degen = mclGetUninitializedArray();
    mxArray * del = mclGetUninitializedArray();
    mxArray * delbnd = mclGetUninitializedArray();
    mxArray * delta = mclGetUninitializedArray();
    mxArray * diff = mclGetUninitializedArray();
    mxArray * dnewt = mclGetUninitializedArray();
    mxArray * dv = mclGetUninitializedArray();
    mxArray * errmsg = mclGetUninitializedArray();
    mxArray * ex = mclGetUninitializedArray();
    mxArray * fdata = mclGetUninitializedArray();
    mxArray * figtr = mclGetUninitializedArray();
    mxArray * findiffevals = mclGetUninitializedArray();
    mxArray * formatstr = mclGetUninitializedArray();
    mxArray * g = mclGetUninitializedArray();
    mxArray * gopt = mclGetUninitializedArray();
    mxArray * grad = mclGetUninitializedArray();
    mxArray * gradflag = mclGetUninitializedArray();
    mxArray * group = mclGetUninitializedArray();
    mxArray * header = mclGetUninitializedArray();
    mxArray * it = mclGetUninitializedArray();
    mxArray * itb = mclGetUninitializedArray();
    mxArray * kmax = mclGetUninitializedArray();
    mxArray * lsotframe = mclGetUninitializedArray();
    mxArray * maxfunevals = mclGetUninitializedArray();
    mxArray * mm = mclGetUninitializedArray();
    mxArray * mtxmpy = mclGetUninitializedArray();
    mxArray * n = mclGetUninitializedArray();
    mxArray * nbnds = mclGetUninitializedArray();
    mxArray * newA = mclGetUninitializedArray();
    mxArray * newfvec = mclGetUninitializedArray();
    mxArray * newgrad = mclGetUninitializedArray();
    mxArray * newval = mclGetUninitializedArray();
    mxArray * newx = mclGetUninitializedArray();
    mxArray * npcg = mclGetUninitializedArray();
    mxArray * nrmsx = mclGetUninitializedArray();
    mxArray * numFunEvals = mclGetUninitializedArray();
    mxArray * numGradEvals = mclGetUninitializedArray();
    mxArray * numberOfVariables = mclGetUninitializedArray();
    mxArray * oposdef = mclGetUninitializedArray();
    mxArray * optnrm = mclGetUninitializedArray();
    mxArray * oval = mclGetUninitializedArray();
    mxArray * p = mclGetUninitializedArray();
    mxArray * pcflags = mclGetUninitializedArray();
    mxArray * pcgit = mclGetUninitializedArray();
    mxArray * pcgtol = mclGetUninitializedArray();
    mxArray * pcmtx = mclGetUninitializedArray();
    mxArray * pert = mclGetUninitializedArray();
    mxArray * posdef = mclGetUninitializedArray();
    mxArray * pp = mclGetUninitializedArray();
    mxArray * prev_diff = mclGetUninitializedArray();
    mxArray * qp = mclGetUninitializedArray();
    mxArray * r = mclGetUninitializedArray();
    mxArray * ratio = mclGetUninitializedArray();
    mxArray * showstat = mclGetUninitializedArray();
    mxArray * snod = mclGetUninitializedArray();
    mxArray * sx = mclGetUninitializedArray();
    mxArray * theta = mclGetUninitializedArray();
    mxArray * tol = mclGetUninitializedArray();
    mxArray * tol1 = mclGetUninitializedArray();
    mxArray * tol2 = mclGetUninitializedArray();
    mxArray * totls = mclGetUninitializedArray();
    mxArray * typx = mclGetUninitializedArray();
    mxArray * v = mclGetUninitializedArray();
    mxArray * val = mclGetUninitializedArray();
    mxArray * vdeg = mclGetUninitializedArray();
    mxArray * vflops = mclGetUninitializedArray();
    mxArray * voptnrm = mclGetUninitializedArray();
    mxArray * vpcg = mclGetUninitializedArray();
    mxArray * vpos = mclGetUninitializedArray();
    mxArray * vval = mclGetUninitializedArray();
    mxArray * warnstr = mclGetUninitializedArray();
    mxArray * x = mclGetUninitializedArray();
    mxArray * xold = mclGetUninitializedArray();
    mclValidateInputs(
      "optim/private/snls",
      12,
      &funfcn,
      &xstart,
      &l,
      &u,
      &verb,
      &options,
      &fval,
      &JACval,
      &YDATA,
      &caller,
      &Jstr,
      &computeLambda);
    mclCopyArray(&xstart);
    mclCopyArray(&l);
    mclCopyArray(&u);
    /*
     * %SNLS  Sparse nonlinear least squares solver.
     * %   
     * %   Locate a local solution
     * %   to the box-constrained nonlinear least-squares problem:
     * %
     * %              min { ||F(x)||^2 :  l <= x <= u}.
     * %
     * %   where F:R^n -> R^m, m > n, and || || is the 2-norm.
     * %
     * % x=SNLS(fname,xstart) solves the unconstrained nonlinear least
     * % squares problem. The vector function is 'fname' and the
     * % starting point is xstart. 
     * %
     * % x=SNLS(fname,xstart,options) solves the unconstrained or
     * % box constrained nonlinear least squares problem. Bounds and
     * % parameter settings are handled through the named parameter list
     * % options.
     * %
     * % x=SNLS(fname,xstart,options,Jstr) indicates the structure
     * % of the sparse Jacobian matrix -- sparse finite differencing
     * % is used to compute the Jacobian when Jstr is not empty.
     * %
     * % x=SNLS(fname,xstart,options,Jstr,fdata) supplies a matrix
     * % fdata to the objective function fdata.o%
     * %
     * % [x,fvec] =SNLS(fname,xstart, ...)  returns the final value of the
     * % vector function F.
     * %
     * % [x,fvec,gopt] = SNLS(fname,xstart, ...) returns the first-order
     * % optimality vector.
     * %
     * % [x,fvec,gopt,it] = SNLS(fname,xstart, ...) returns the number of
     * % major iterations used.
     * %
     * % [x,fvec,gopt,it,npcg] =  SNLS(fname,xstart, ...) returns the
     * % total number of CG iterations used.
     * %
     * % [x,fvec,gopt,it,npcg,ex] = SNLS(fname,xstart, ...) returns the 
     * % termination code.
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.16 $  $Date: 1998/09/11 19:29:49 $
     * 
     * %   Extract input parameters etc.
     * l=l(:); u = u(:);
     */
    mlfAssign(&l, mlfIndexRef(l, "(?)", mlfCreateColonIndex()));
    mlfAssign(&u, mlfIndexRef(u, "(?)", mlfCreateColonIndex()));
    /*
     * xcurr=xstart;  % save shape of xstart for user function
     */
    mlfAssign(&xcurr, xstart);
    /*
     * xstart = xstart(:);  % make it a vector
     */
    mlfAssign(&xstart, mlfIndexRef(xstart, "(?)", mlfCreateColonIndex()));
    /*
     * 
     * %   Initialize
     * dnewt = []; fdata = [];
     */
    mlfAssign(&dnewt, mclCreateEmptyArray());
    mlfAssign(&fdata, mclCreateEmptyArray());
    /*
     * A = []; 
     */
    mlfAssign(&A, mclCreateEmptyArray());
    /*
     * fvec = []; 
     */
    mlfAssign(fvec, mclCreateEmptyArray());
    /*
     * gopt = [];
     */
    mlfAssign(&gopt, mclCreateEmptyArray());
    /*
     * n = length(xstart); 
     */
    mlfAssign(&n, mlfLength(xstart));
    /*
     * it= 1;  
     */
    mlfAssign(&it, mlfScalar(1.0));
    /*
     * numFunEvals = 1;  % done in calling function lsqnonlin
     */
    mlfAssign(&numFunEvals, mlfScalar(1.0));
    /*
     * numGradEvals = 1; % done in calling function
     */
    mlfAssign(&numGradEvals, mlfScalar(1.0));
    /*
     * totls = 0;
     */
    mlfAssign(&totls, mlfScalar(0.0));
    /*
     * header = sprintf(['\n                                         Norm of      First-order \n',...
     */
    mlfAssign(
      &header,
      mlfSprintf(
        NULL,
        mlfHorzcat(
          mxCreateString(
            "\\n                                  "
            "       Norm of      First-order \\n"),
          mxCreateString(
            " Iteration  Func-count     f(x)          s"
            "tep          optimality   CG-iterations"),
          NULL),
        NULL));
    /*
     * ' Iteration  Func-count     f(x)          step          optimality   CG-iterations']);
     * formatstr = ' %5.0f      %5.0f   %13.6g  %13.6g   %12.3g      %7.0f';
     */
    mlfAssign(
      &formatstr,
      mxCreateString(" %5.0f      %5.0f   %13.6g  %13.6g   %12.3g      %7.0f"));
    /*
     * 
     * if n == 0, 
     */
    if (mlfTobool(mlfEq(n, mlfScalar(0.0)))) {
        /*
         * error('n must be positive'), 
         */
        mlfError(mxCreateString("n must be positive"));
    /*
     * end
     */
    }
    /*
     * 
     * if isempty(l), 
     */
    if (mlfTobool(mlfIsempty(l))) {
        /*
         * l = -inf*ones(n,1); 
         */
        mlfAssign(
          &l, mlfMtimes(mlfUminus(mlfInf()), mlfOnes(n, mlfScalar(1.0), NULL)));
    /*
     * end
     */
    }
    /*
     * if isempty(u), 
     */
    if (mlfTobool(mlfIsempty(u))) {
        /*
         * u = inf*ones(n,1); 
         */
        mlfAssign(&u, mlfMtimes(mlfInf(), mlfOnes(n, mlfScalar(1.0), NULL)));
    /*
     * end
     */
    }
    /*
     * arg = (u >= 1e10); 
     */
    mlfAssign(&arg, mlfGe(u, mlfScalar(1e10)));
    /*
     * arg2 = (l <= -1e10);
     */
    mlfAssign(&arg2, mlfLe(l, mlfScalar(-1e10)));
    /*
     * u(arg) = inf;
     */
    mlfIndexAssign(&u, "(?)", arg, mlfInf());
    /*
     * l(arg2) = -inf;
     */
    mlfIndexAssign(&l, "(?)", arg2, mlfUminus(mlfInf()));
    /*
     * if any(u == l) 
     */
    if (mlfTobool(mlfAny(mlfEq(u, l), NULL))) {
        /*
         * errmsg=sprintf('%s',...
         */
        mlfAssign(
          &errmsg,
          mlfSprintf(
            NULL,
            mxCreateString("%s"),
            mxCreateString("Equal upper and lower bounds not permitted."),
            NULL));
        /*
         * 'Equal upper and lower bounds not permitted.');
         * error(errmsg)
         */
        mlfError(errmsg);
    /*
     * elseif min(u-l) <= 0
     */
    } else if (mlfTobool(
                 mlfLe(
                   mlfMin(NULL, mlfMinus(u, l), NULL, NULL), mlfScalar(0.0)))) {
        /*
         * error('Inconsistent bounds.')
         */
        mlfError(mxCreateString("Inconsistent bounds."));
    /*
     * end
     */
    }
    /*
     * if min(min(u-xstart),min(xstart-l)) < 0
     */
    if (mlfTobool(
          mlfLt(
            mlfMin(
              NULL,
              mlfMin(NULL, mlfMinus(u, xstart), NULL, NULL),
              mlfMin(NULL, mlfMinus(xstart, l), NULL, NULL),
              NULL),
            mlfScalar(0.0)))) {
        /*
         * xstart = startx(u,l);    
         */
        mlfAssign(&xstart, mlfOptim_private_startx(u, l));
    /*
     * end
     */
    }
    /*
     * 
     * %
     * % get options out
     * active_tol = optimget(options,'ActiveConstrTol',sqrt(eps));
     */
    mlfAssign(
      &active_tol,
      mlfOptimget(
        options, mxCreateString("ActiveConstrTol"), mlfSqrt(mlfEps())));
    /*
     * gradflag =  strcmp(optimget(options,'Jacobian'),'on');
     */
    mlfAssign(
      &gradflag,
      mlfStrcmp(
        mlfOptimget(options, mxCreateString("Jacobian"), NULL),
        mxCreateString("on")));
    /*
     * typx = optimget(options,'typicalx') ;
     */
    mlfAssign(&typx, mlfOptimget(options, mxCreateString("typicalx"), NULL));
    /*
     * % In case the defaults were gathered from calling: optimset('quadprog'):
     * numberOfVariables = n;
     */
    mlfAssign(&numberOfVariables, n);
    /*