sfdnls.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 "sfdnls.h"

/*
 * The function "Msfdnls" is the implementation version of the "sfdnls"
 * M-function from file "C:\MATLABR11\toolbox\optim\sfdnls.m" (lines 1-89). 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[J, ncol] = sfdnls(xcurr,valx,Jstr,group,alpha,fun,YDATA,varargin)
 */
static mxArray * Msfdnls(mxArray * * ncol,
                         int nargout_,
                         mxArray * xcurr,
                         mxArray * valx,
                         mxArray * Jstr,
                         mxArray * group,
                         mxArray * alpha,
                         mxArray * fun,
                         mxArray * YDATA,
                         mxArray * varargin) {
    mxArray * J = mclGetUninitializedArray();
    mxArray * A = mclGetUninitializedArray();
    mxArray * cols = mclGetUninitializedArray();
    mxArray * d = mclGetUninitializedArray();
    mxArray * epsi = mclGetUninitializedArray();
    mxArray * i = mclGetUninitializedArray();
    mxArray * ind = mclGetUninitializedArray();
    mclForLoopIterator iterator_0;
    mxArray * j = mclGetUninitializedArray();
    mxArray * k = mclGetUninitializedArray();
    mxArray * lpoint = mclGetUninitializedArray();
    mxArray * m = mclGetUninitializedArray();
    mxArray * n = mclGetUninitializedArray();
    mxArray * nargin_ = mclGetUninitializedArray();
    mxArray * p = mclGetUninitializedArray();
    mxArray * v = mclGetUninitializedArray();
    mxArray * val = mclGetUninitializedArray();
    mxArray * w = mclGetUninitializedArray();
    mxArray * x = mclGetUninitializedArray();
    mxArray * xnrm = mclGetUninitializedArray();
    mxArray * y = mclGetUninitializedArray();
    mlfAssign(
      &nargin_,
      mlfNargin(
        1, xcurr, valx, Jstr, group, alpha, fun, YDATA, varargin, NULL));
    mclValidateInputs(
      "sfdnls", 7, &xcurr, &valx, &Jstr, &group, &alpha, &fun, &YDATA);
    mclCopyArray(&xcurr);
    mclCopyArray(&alpha);
    mclCopyArray(&YDATA);
    /*
     * %SFDNLS    Sparse Jacobian via finite differences
     * %
     * %
     * % J = sfdnls(x,valx,J,group,[],fun) returns the
     * % sparse finite difference approximation J of a Jacobian matrix
     * % of function 'fun'  at current point x.
     * % Vector group indicates how to use sparse finite differencing:
     * % group(i) = j means that column i belongs to group (or color) j.
     * % Each group (or color) corresponds to a function difference.
     * % varargin are extra parameters (possibly) needed by function 'fun'.
     * %
     * % J = sfdnls(x,valx,J,group,fdata,fun,alpha) overrides the default
     * % finite differencing stepsize.
     * %
     * %[J, ncol] = sfdnls(...) returns the number of function evaluations used
     * % in ncol.
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.6 $  $Date: 1998/07/22 23:08:33 $
     * 
     * %
     * if nargin < 6
     */
    if (mlfTobool(mlfLt(nargin_, mlfScalar(6.0)))) {
        /*
         * error('SFDNLS requires six arguments')
         */
        mlfError(mxCreateString("SFDNLS requires six arguments"));
    /*
     * elseif nargin < 7
     */
    } else if (mlfTobool(mlfLt(nargin_, mlfScalar(7.0)))) {
        /*
         * YDATA = [];
         */
        mlfAssign(&YDATA, mclCreateEmptyArray());
    /*
     * end
     */
    }
    /*
     * 
     * x = xcurr(:); % make it a vector
     */
    mlfAssign(&x, mlfIndexRef(xcurr, "(?)", mlfCreateColonIndex()));
    /*
     * [m,n] = size(Jstr); 
     */
    mlfSize(mlfVarargout(&m, &n, NULL), Jstr, NULL);
    /*
     * ncol = max(group); epsi = sqrt(eps);
     */
    mlfAssign(ncol, mlfMax(NULL, group, NULL, NULL));
    mlfAssign(&epsi, mlfSqrt(mlfEps()));
    /*
     * if isempty(alpha), 
     */
    if (mlfTobool(mlfIsempty(alpha))) {
        /*
         * alpha = ones(ncol,1)*sqrt(eps); 
         */
        mlfAssign(
          &alpha,
          mlfMtimes(mlfOnes(*ncol, mlfScalar(1.0), NULL), mlfSqrt(mlfEps())));
    /*
     * end
     */
    }
    /*
     * J = spones(Jstr); d = zeros(n,1);
     */
    mlfAssign(&J, mlfSpones(Jstr));
    mlfAssign(&d, mlfZeros(n, mlfScalar(1.0), NULL));
    /*
     * if ncol < n
     */
    if (mlfTobool(mlfLt(*ncol, n))) {
        /*
         * for k = 1:ncol
         */
        for (mclForStart(&iterator_0, mlfScalar(1.0), *ncol, NULL);
             mclForNext(&iterator_0, &k);
             ) {
            /*
             * d = (group == k);
             */
            mlfAssign(&d, mlfEq(group, k));
            /*
             * if nargin < 5
             */
            if (mlfTobool(mlfLt(nargin_, mlfScalar(5.0)))) {
                /*
                 * xnrm = norm(x(d));
                 */
                mlfAssign(&xnrm, mlfNorm(mlfIndexRef(x, "(?)", d), NULL));
                /*
                 * xnrm = max(xnrm,1);
                 */
                mlfAssign(&xnrm, mlfMax(NULL, xnrm, mlfScalar(1.0), NULL));
                /*
                 * alpha(k) = alpha(k)*xnrm;
                 */
                mlfIndexAssign(
                  &alpha,
                  "(?)",
                  k,
                  mlfMtimes(mlfIndexRef(alpha, "(?)", k), xnrm));
            /*
             * end
             */
            }
            /*
             * y = x + alpha(k)*d;
             */
            mlfAssign(
              &y, mlfPlus(x, mlfMtimes(mlfIndexRef(alpha, "(?)", k), d)));
            /*
             * 
             * %v = feval(fun,y,fdata);
             * xcurr(:) = y;  % reshape for userfunction
             */
            mlfIndexAssign(&xcurr, "(?)", mlfCreateColonIndex(), y);
            /*
             * v = feval(fun,xcurr,varargin{:});
             */
            mlfAssign(
              &v,
              mlfFeval(
                mclValueVarargout(),
                mclFevalLookup(fun, 0, NULL),
                xcurr,
                mlfIndexRef(varargin, "{?}", mlfCreateColonIndex()),
                NULL));
            /*
             * if ~isempty(YDATA)
             */
            if (mlfTobool(mlfNot(mlfIsempty(YDATA)))) {
                /*
                 * v = v - YDATA;
                 */
                mlfAssign(&v, mlfMinus(v, YDATA));
            /*
             * end
             */
            }
            /*
             * v = v(:);
             */
            mlfAssign(&v, mlfIndexRef(v, "(?)", mlfCreateColonIndex()));
            /*
             * 
             * w = (v-valx)/alpha(k);
             */
            mlfAssign(
              &w, mlfMrdivide(mlfMinus(v, valx), mlfIndexRef(alpha, "(?)", k)));
            /*
             * cols = find(d); 
             */
            mlfAssign(&cols, mlfFind(NULL, NULL, d));
            /*
             * 
             * lpoint = length(cols);
             */
            mlfAssign(&lpoint, mlfLength(cols));
            /*
             * A = sparse(m,n);
             */
            mlfAssign(&A, mlfSparse(m, n, NULL, NULL, NULL, NULL));
            /*
             * A(:,cols) = J(:,cols);
             */
            mlfIndexAssign(
              &A,
              "(?,?)",
              mlfCreateColonIndex(),
              cols,
              mlfIndexRef(J, "(?,?)", mlfCreateColonIndex(), cols));
            /*
             * J(:,cols) = J(:,cols) - A(:,cols);
             */
            mlfIndexAssign(
              &J,
              "(?,?)",
              mlfCreateColonIndex(),
              cols,
              mlfMinus(
                mlfIndexRef(J, "(?,?)", mlfCreateColonIndex(), cols),
                mlfIndexRef(A, "(?,?)", mlfCreateColonIndex(), cols)));
            /*
             * [i,j,val] = find(A);
             */
            mlfAssign(&i, mlfFind(&j, &val, A));
            /*
             * [p,ind] = sort(i);
             */
            mlfAssign(&p, mlfSort(&ind, i, NULL));
            /*
             * val(ind) = w(p);
             */
            mlfIndexAssign(&val, "(?)", ind, mlfIndexRef(w, "(?)", p));