optim_private_pcgr.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_pcgr.h"

static mxArray * Moptim_private_pcgr(mxArray * * posdef,
                                     mxArray * * k,
                                     int nargout_,
                                     mxArray * DM,
                                     mxArray * DG,
                                     mxArray * g,
                                     mxArray * kmax,
                                     mxArray * tol,
                                     mxArray * mtxmpy,
                                     mxArray * fdata,
                                     mxArray * H,
                                     mxArray * R,
                                     mxArray * pR);
static mxArray * Mpcgr_preproj(int nargout_,
                               mxArray * r,
                               mxArray * RPCMTX,
                               mxArray * ppvec);
static mxArray * mlfPcgr_preproj(mxArray * r,
                                 mxArray * RPCMTX,
                                 mxArray * ppvec);
static void mlxPcgr_preproj(int nlhs,
                            mxArray * plhs[],
                            int nrhs,
                            mxArray * prhs[]);

static mlfFunctionTableEntry local_function_table_[1]
  = { { "preproj", mlxPcgr_preproj, 3, 1 } };

/*
 * The function "Moptim_private_pcgr" is the implementation version of the
 * "optim/private/pcgr" M-function from file
 * "C:\MATLABR11\toolbox\optim\private\pcgr.m" (lines 1-74). 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[p,posdef,k] = pcgr(DM,DG,g,kmax,tol,mtxmpy,fdata,H,R,pR);
 */
static mxArray * Moptim_private_pcgr(mxArray * * posdef,
                                     mxArray * * k,
                                     int nargout_,
                                     mxArray * DM,
                                     mxArray * DG,
                                     mxArray * g,
                                     mxArray * kmax,
                                     mxArray * tol,
                                     mxArray * mtxmpy,
                                     mxArray * fdata,
                                     mxArray * H,
                                     mxArray * R,
                                     mxArray * pR) {
    mxArray * p = mclGetUninitializedArray();
    mxArray * alpha = mclGetUninitializedArray();
    mxArray * beta = mclGetUninitializedArray();
    mxArray * d = mclGetUninitializedArray();
    mxArray * denom = mclGetUninitializedArray();
    mxArray * inner1 = mclGetUninitializedArray();
    mxArray * inner2 = mclGetUninitializedArray();
    mclForLoopIterator iterator_0;
    mxArray * m = mclGetUninitializedArray();
    mxArray * n = mclGetUninitializedArray();
    mxArray * r = mclGetUninitializedArray();
    mxArray * stoptol = mclGetUninitializedArray();
    mxArray * val = mclGetUninitializedArray();
    mxArray * w = mclGetUninitializedArray();
    mxArray * ww = mclGetUninitializedArray();
    mxArray * z = mclGetUninitializedArray();
    mxArray * znrm = mclGetUninitializedArray();
    mclValidateInputs(
      "optim/private/pcgr",
      10,
      &DM,
      &DG,
      &g,
      &kmax,
      &tol,
      &mtxmpy,
      &fdata,
      &H,
      &R,
      &pR);
    mclCopyArray(&kmax);
    /*
     * %PCGR	Preconditioned conjugate gradients
     * %
     * % [p,posdef,k] = PCGR(DM,DG,g,kmax,tol,mtxmpy,fdata,H,R,pR) apply
     * % a preconditioned conjugate gradient procedure to the quadratic
     * % 
     * %         q(p) = .5p'Mp + g'p, where
     * %
     * % M = DM*H*DM + DG. kmax is a bound on the number of permitted
     * % CG-iterations, tol is a stopping tolerance on the residual (default
     * % is tol = .1), mtxmpy is the function that computes products
     * % with the Hessian matrix H (possible using data fdata),  
     * % and R is the cholesky factor of the preconditioner (transpose) of
     * % M. So, R'R approximates M(pR,pR), where pR is a permutation vector.
     * % On ouput p is the computed direction, posdef = 1 implies
     * % only positive curvature (in M) has been detected; posdef = 0
     * % implies p is a direction of negative curvature (for M).
     * % Output parameter k is the number of CG-iterations used (which
     * % corresponds to the number of multiplications with H).
     * %
     * 
     * %   Copyright (c) 1990-98 by The MathWorks, Inc.
     * %   $Revision: 1.3 $  $Date: 1998/08/17 19:42:07 $
     * 
     * % Initializations.
     * n = length(DG);
     */
    mlfAssign(&n, mlfLength(DG));
    /*
     * r = -g; 
     */
    mlfAssign(&r, mlfUminus(g));
    /*
     * p = zeros(n,1); 
     */
    mlfAssign(&p, mlfZeros(n, mlfScalar(1.0), NULL));
    /*
     * val = 0;
     */
    mlfAssign(&val, mlfScalar(0.0));
    /*
     * m = 0;
     */
    mlfAssign(&m, mlfScalar(0.0));
    /*
     * 
     * % Precondition .
     * z = preproj(r,R,pR);
     */
    mlfAssign(&z, mlfPcgr_preproj(r, R, pR));
    /*
     * znrm = norm(z); 
     */
    mlfAssign(&znrm, mlfNorm(z, NULL));
    /*
     * stoptol = tol*znrm;
     */
    mlfAssign(&stoptol, mlfMtimes(tol, znrm));
    /*
     * inner2 = 0; 
     */
    mlfAssign(&inner2, mlfScalar(0.0));
    /*
     * inner1 = r'*z;
     */
    mlfAssign(&inner1, mlfMtimes(mlfCtranspose(r), z));
    /*
     * posdef = 1;
     */
    mlfAssign(posdef, mlfScalar(1.0));
    /*
     * 
     * kmax = max(kmax,1);  % kmax must be at least 1
     */
    mlfAssign(&kmax, mlfMax(NULL, kmax, mlfScalar(1.0), NULL));
    /*
     * % PRIMARY LOOP.
     * for k = 1:kmax
     */
    for (mclForStart(&iterator_0, mlfScalar(1.0), kmax, NULL);
         mclForNext(&iterator_0, k);
         ) {
        /*
         * if k==1
         */
        if (mlfTobool(mlfEq(*k, mlfScalar(1.0)))) {
            /*
             * d = z;
             */
            mlfAssign(&d, z);
        /*
         * else
         */
        } else {
            /*
             * beta = inner1/inner2;
             */
            mlfAssign(&beta, mlfMrdivide(inner1, inner2));
            /*
             * d = z + beta*d;
             */
            mlfAssign(&d, mlfPlus(z, mlfMtimes(beta, d)));
        /*
         * end
         */
        }
        /*
         * ww = DM*d; 
         */
        mlfAssign(&ww, mlfMtimes(DM, d));
        /*
         * w = feval(mtxmpy,ww,H,fdata);
         */
        mlfAssign(
          &w,
          mlfFeval(
            mclValueVarargout(),
            mclFevalLookup(mtxmpy, 1, local_function_table_),
            ww,
            H,
            fdata,
            NULL));
        /*
         * ww = DM*w +DG*d;
         */
        mlfAssign(&ww, mlfPlus(mlfMtimes(DM, w), mlfMtimes(DG, d)));
        /*
         * denom = d'*ww;
         */
        mlfAssign(&denom, mlfMtimes(mlfCtranspose(d), ww));
        /*
         * if denom <= 0
         */
        if (mlfTobool(mlfLe(denom, mlfScalar(0.0)))) {
            /*
             * p = d/norm(d);
             */
            mlfAssign(&p, mlfMrdivide(d, mlfNorm(d, NULL)));
            /*
             * posdef = 0;
             */
            mlfAssign(posdef, mlfScalar(0.0));
            /*
             * break
             */
            mclDestroyForLoopIterator(&iterator_0);
            break;
        /*
         * else
         */
        } else {
            /*
             * alpha = inner1/denom;
             */
            mlfAssign(&alpha, mlfMrdivide(inner1, denom));
            /*
             * p = p+ alpha*d;
             */
            mlfAssign(&p, mlfPlus(p, mlfMtimes(alpha, d)));
            /*
             * r = r - alpha*ww;
             */
            mlfAssign(&r, mlfMinus(r, mlfMtimes(alpha, ww)));
        /*
         * end
         */
        }
        /*
         * z = preproj(r,R,pR);
         */
        mlfAssign(&z, mlfPcgr_preproj(r, R, pR));
        /*
         * 
         * % Exit?
         * if norm(z) <= stoptol 
         */
        if (mlfTobool(mlfLe(mlfNorm(z, NULL), stoptol))) {
            /*
             * break; 
             */
            mclDestroyForLoopIterator(&iterator_0);
            break;
        /*
         * end
         */
        }
        /*
         * inner2 = inner1;
         */
        mlfAssign(&inner2, inner1);
        /*
         * inner1 = r'*z;
         */
        mlfAssign(&inner1, mlfMtimes(mlfCtranspose(r), z));
    /*
     * end
     */
    }
    mclValidateOutputs("optim/private/pcgr", 3, nargout_, &p, posdef, k);
    mxDestroyArray(alpha);
    mxDestroyArray(beta);
    mxDestroyArray(d);
    mxDestroyArray(denom);
    mxDestroyArray(inner1);
    mxDestroyArray(inner2);
    mxDestroyArray(kmax);
    mxDestroyArray(m);
    mxDestroyArray(n);
    mxDestroyArray(r);
    mxDestroyArray(stoptol);
    mxDestroyArray(val);
    mxDestroyArray(w);
    mxDestroyArray(ww);
    mxDestroyArray(z);
    mxDestroyArray(znrm);
    /*
     * 
     * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     * 
     */
    return p;
}

/*
 * The function "mlfOptim_private_pcgr" contains the normal interface for the
 * "optim/private/pcgr" M-function from file