eqsolve.src

来自「没有说明」· SRC 代码 · 共 549 行 · 第 1/2 页

        end;
    endif;

/*
** Initialize defaults.
**
*/
    retcode = 0;

    eta = sqrt(__macheps);
    eqmtol = __macheps^(2/3);

    if eqstol == 0;
        eqstol = __macheps^(2/3);
    endif;

    if eqtol == 0;
        eqtol = __macheps^(1/3);
    endif;


    if eqtypx == 0;
        sx = ones(n,1);
        typx = sx;
    elseif rows(sx) == n;
        typx = missrv(miss(eqtypx,0),1);
        sx = 1/typx;
    endif;

    if eqtypf == 0;
        sf = ones(n,1);
        typf = sf;
    elseif rows(sf) == n;
        typf = missrv(miss(eqtypf,0),1);
        sf = 1/typf;
    endif;

    sfsq = sf^2;

    maxstep = (10^3)*maxc(sqrt(moment(sx.*x0,0))'|sqrt(moment(sx,0))');

    if eqaltnam $== 0;
        varlist = zeros(n,1) $+ "X" $+ ftocv(seqa(1,1,n), eqvpad*(floor(log(n)
            )+1),0);
    else;
        if rows(eqaltnam) == n;
            varlist = eqaltnam;
        elseif cols(eqaltnam) == n;
            varlist = eqaltnam';
        else;
            errorlog "\nERROR:  Wrong number of rows in __altnam.\n";
            end;
        endif;
    endif;
    eqlist = 0 $+ "F" $+ ftocv(seqa(1,1,n), eqvpad*(floor(log(n))+1),0) $+
        "(X)";
    fmtstr = { "-*.*lf" 15 8, "*.*lf" 18 8, "*.*lf" 22 8, "*.*lf" 22 8 };

    if eqjac;
        jacob = eqjac;      /* Assign pointer to jacob */
        local jacob:proc;           /* Assume jacob is procedure pointer  */
        jc = jacob(x0);
    else;
        jc = _eqs_jacob(&f,x0,fvc,0,typx,eta);
    endif;

    gc = jc'fvc;
    xp = x0;
    fc = 0.5*moment(sf.*fvc,0);
    retcode = 0;

    consecmx = 0;
    itncount = 0;

    if eqiterinfo;
        print chrs(45*ones(80,1));
        print "              ROOTS";;
        print chrs(32*ones(32,1)) "      F(ROOTS)";
        print chrs(45*ones(80,1));
        print;
    endif;

    do while retcode == 0;

        p = -fvc/jc;
        gc = jc'fvc;
        xc = xp;

        local alpha,newtlen,rlength,minlam,lambda,ab,inslope,
            ltemp, one,two,three,fpprev,a,b,disc,lprev, notdone,
            tcode;

    /* Initialize Parameters */

        maxtaken = 0;
        alpha = 1e-4;

        newtlen = sqrt(moment(sx.*p,0));
        if newtlen > maxstep;
            p = p*(maxstep/newtlen);
            newtlen = maxstep;
        endif;

        inslope = gc'p;
        rlength = maxc(abs(p) ./ (maxc(xc'|(typx')) ));

        if rlength == 0;
            xp = xc;
            fp = fc;
        else;
            minlam = eqstol/rlength;

            lambda = 1;
            notdone = 1;

            do while notdone;
                xp = xc + lambda*p;
                fp = 0.5*moment(sf.*f(xp),0);       /* objective function  */

                if fp <= fc + alpha*lambda*inslope;
                    if lambda == 1 and newtlen > .99*maxstep;
                       maxtaken = 1;
                    endif;
                    notdone = 0;
                elseif lambda < minlam;
                    xp = xc;
                    notdone = 0;
                else;
                    if lambda == 1;
                            /* Quadratic step */
                        ltemp = -inslope/(2*(fp-fc-inslope));
                    else;
                        one = 1/(lambda - lprev);
                        two = ((1/lambda^2 ~ -1/lprev^2) |
                            (-lprev/lambda^2 ~ lambda/lprev^2));
                        three = fp - fc - lambda *inslope | fpprev - fc
                            -lprev*inslope;

                        ab = one*two*three;
                        a = ab[1,1];
                        b = ab[2,1];
                        disc = b^2 - 3*a*inslope;

                        if ndpchk(8);
                            break;

                        elseif a == 0;
                            ltemp = -inslope/(2*b);
                        elseif disc <= 0 or ndpchk(8);
                            break;
                        else;
                            ltemp = (-b+sqrt(disc))/(3*a);
                        endif;
                        if ltemp > 0.5*lambda;
                            ltemp = 0.5*lambda;
                        endif;
                    endif;
                    lprev = lambda;
                    fpprev = fp;

                    if ltemp <= 0.1*lambda;
                        lambda = 0.1*lambda;
                    else;
                        lambda = ltemp;
                    endif;
                endif;
            endo;

        endif;

        fvp = f(xp);

        if eqjac;
            jc = jacob(xp);
        else;
            jc = _eqs_jacob(&f,xp,fvp,0,typx,eta);
        endif;

/*
** Check to see if any stopping criteria have been met.
*/

        if _eqs_norm2(fvp.*sf) <= eqtol;
            retcode = 1;    /* A proper termination */
        elseif abs(xp-xc)./maxc(abs(xp)'|typx') <= eqstol;
            retcode = 2;    /* Algorithm may have bogged down */
        elseif retcode == 1;
            retcode = 3;
        elseif itncount >= eqmaxit;
            retcode = 4;
        elseif maxtaken == 1;
            consecmx = consecmx+1;
            if consecmx == 20;
                retcode = 5;
            endif;
        else;
            consecmx = 0;
            gp = jc'.*sfsq*fvp;
            if maxc(abs((gp.* maxc(xp'|typx'))/maxc(fp|1/sf))) <= eqmtol;
                retcode = 6;
            endif;
        endif;
        itncount = itncount + 1;

        if eqiterinfo > 0;
            print "Iteration #" ftos(itncount,"%*.*lf",1,0);
            call printfm(varlist~xp~eqlist~fvp,0~1~0~1,fmtstr);
            print;
            print;
        endif;

        fc = fp;
        fvc = fvp;

    endo;


retp(xp,retcode);
endp;

proc _eqs_norm2(x);
    retp(0.5*sqrt(x'x));
endp;




proc _eqs_jacob(f,x,f0,dh,typx,eta);

    local rows_f,rows_x,jacobian,i,xc,xp;
    local f:proc;

    rows_f = rows(f0);
    rows_x = rows(x);
    jacobian = zeros(rows_f,rows_x);

    /* Computation of stepsize (dh) for gradient if dh == 0 */

    if dh == 0;
        dh = eta*maxc(abs(x')|typx');
        dh = dh + (dh == 0);
    endif;

    xc = x + dh;
    dh = xc - x;        /* This increases precision slightly */
    xp = x + eye(rows_x).*dh;

    /* Calculate forward difference */
    i = 1;
    do until i > rows_x;
        jacobian[.,i] = (f(xp[.,i]) - f0);
        i = i + 1;
    endo;
    retp(jacobian./(dh'));
endp;




proc (0) = eqSolveSet;

    gausset;

    _eqs_typicalX = 0;
    _eqs_typicalF = 0;
    _eqs_IterInfo = 0;
    _eqs_JacobianProc = 0;
    _eqs_MaxIters = 100;
    _eqs_StepTol = __macheps^(2/3);

endp;