asa_cg.cpp

数学计算程序

            if ( phi <= fquad )
            {
                if ( Parm->PrintLevel >= 4 )
                {
                    printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n",
                            alpha, phi, fquad) ;
                }
                if ( asa_Wolfe (alpha, phi, dphi, Com) )
                {
                    status = 0 ;
                    goto Exit ;
                }
            }
        }
        if ( phi <= Com->fpert )
        {
            a = alpha ;
            dphia = dphi ;
        }
        else
        {
            /* contraction phase, only break at termination or Secant step */
            b = alpha ;
            while ( TRUE )
            {
                alpha = .5*(a+b) ;
                nshrink++ ;
                if ( nshrink > Parm->nexpand )
                {
                    status = 6 ;
                    goto Exit ;
                }
                asa_step (xtemp, x, d, alpha, nfree) ;
                asa_g (gtemp, xtemp, Com) ;
                dphi = asa_dot (gtemp, d, nfree) ;
                if ( dphi >= ZERO ) goto Secant ;
                phi = asa_f (xtemp, Com) ;
                if ( Parm->PrintLevel >= 4 )
                {
                    printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: "
                            "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ;
                }
                if ( Com->QuadOK && (phi <= fquad) )
                {
                    if ( asa_Wolfe (alpha, phi, dphi, Com) )
                    {
                        status = 0 ;
                        goto Exit ;
                    }
                }
                if ( phi <= Com->fpert )
                {
                    a = alpha ;
                    dphia = dphi ;
                }
                else
                {
                    b = alpha ;
                }
            }
        }

/* expansion phase */

        ngrow++ ;
        if ( ngrow > Parm->nexpand )
        {
            status = 3 ;
            goto Exit ;
        }
        alphaold = alpha ;
        alpha = MIN (rho*alpha, minstep) ;
        if ( alpha != alphaold )
        {
            asa_step (xtemp, x, d, alpha, nfree) ;
            asa_g (gtemp, xtemp, Com) ;
            dphi = asa_dot (gtemp, d, nfree) ;
            if ( Parm->PrintLevel >= 4 )
            {
                printf ("expand,   a: %14.6e alpha: %14.6e phi: "
                         "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ;
            }
        }
        else /* a new constraint is active */
        {
            do /* while statement */
            {
                alphaold = alpha ;
                phiold = phi ;
                if ( alpha < maxstep )
                {
                    alpha = rho*alphaold ;
                    asa_project (xtemp, x, d, alpha, Com) ;
                    phi = asa_f (xtemp, Com) ;
                }
            } while ( phi < phiold ) ;
            if ( alphaold == minstep )
            {
                asa_step (xtemp, x, d, minstep, nfree) ;
                status = -3 ;
            }
            else
            {
                asa_project (xtemp, x, d, alphaold, Com) ;
                asa_g (gtemp, xtemp, Com) ;
                status = -4 ;
            }
            phi = phiold ;
            goto Exit ;
        }
    }

Secant:
    b = alpha ;
    dphib = dphi ;
    if ( Com->QuadOK )
    {
        phi = asa_f (xtemp, Com) ;
        if ( ngrow + nshrink == 0 ) fquad = MIN (phi, Com->f0) ;
        if ( phi <= fquad )
        {
            if ( asa_Wolfe (alpha, phi, dphi, Com) )
            {
                status = 0 ;
                goto Exit ;
            }
        }
    }
    nsecant = Parm->nsecant ;
    for (iter = 1; iter <= nsecant; iter++)
    {
        if ( Parm->PrintLevel >= 4 )
        {
            printf ("secant, a: %14.6e b: %14.6e da: %14.6e db: %14.6e\n",
                     a, b, dphia, dphib) ;
        }
        width = Parm->gamma*(b - a) ;
        if ( -dphia <= dphib ) alpha = a - (a-b)*(dphia/(dphia-dphib)) ;
        else                   alpha = b - (a-b)*(dphib/(dphia-dphib)) ;
        c = alpha ;
        a0 = a ;
        b0 = b ;
        da0 = dphia ;
        db0 = dphib ;
        status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi,
                    &dphi, Com) ;
        if ( status >= 0 ) goto Exit ;
        else if ( status == -2 )
        {
            if ( c == a )
            {
                if ( dphi > da0 ) alpha = c - (c-a0)*(dphi/(dphi-da0)) ;
                else              alpha = a ;
            }
            else
            {
                if ( dphi < db0 ) alpha = c - (c-b0)*(dphi/(dphi-db0)) ;
                else              alpha = b ;
            }
            if ( (alpha > a) && (alpha < b) )
            {
                if ( Parm->PrintLevel >= 4 ) printf ("2nd secant\n") ;
                status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi,
                          &dphi, Com) ;
                if ( status >= 0 ) goto Exit ;
            }
        }

/* bisection iteration */

        if ( b-a >= width )
        {
            alpha = .5*(b+a) ;
            if ( Parm->PrintLevel >= 4 ) printf ("bisection\n") ;
            status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi,
                        &dphi, Com) ;
            if ( status >= 0 ) goto Exit ;
        }
        else if ( b <= a )
        {
            status = 7 ;
            goto Exit ;
        }
    }
    status = 4 ;

Exit:
    Com->alpha = alpha ;
    Com->f = phi ;
    Com->df = dphi ;
    return (status) ;
}

/* =========================================================================
   === asa_lineW ===========================================================
   =========================================================================
   Ordinary Wolfe line search routine.
   This routine is identical to asa_line except that the function
   psi [a] = phi [a] - phi [0] - a*delta*dphi [0] is minimized instead of
   the function phi
   ========================================================================= */
int asa_lineW
(
    double       dphi0 , /* function derivative at starting point (alpha = 0) */
    asa_com       *Com   /* cg com structure */
)
{
    int i, iter, nfree, nsecant, nshrink, ngrow, status ;
    double a, dpsia, b, dpsib, c, alpha, phi, dphi, alphaold, phiold,
           a0, da0, b0, db0, width, fquad, rho, psi, dpsi, minstep, maxstep,
           *x, *d, *xtemp, *gtemp ;
    asacg_parm *Parm ;

    nfree = Com->nfree ;
    x = Com->x ;         /* current iterate */
    d = Com->d ;         /* current search direction */
    xtemp = Com->xtemp ; /* x + alpha*d */
    gtemp = Com->gtemp ; /* gradient at x + alpha*d */
    minstep = Com->minstep ;
    maxstep = Com->maxstep ;
    alpha = Com->alpha ;
    if ( alpha > minstep )
    {
        alpha = minstep ;
        Com->QuadOK = FALSE ;
    }
    phi = Com->f ;
    Parm = Com->cgParm ;
    rho = Parm->rho ;
    asa_step (xtemp, x, d, alpha, nfree) ;
    asa_g (gtemp, xtemp, Com) ;
    dphi = asa_dot (gtemp, d, nfree) ;

    /* check if gradient is nan; if so, reduce stepsize */
    if ( dphi != dphi )
    {
        for (i = 0; i < Parm->nexpand; i++)
        {
            alpha /= rho ;
            asa_step (xtemp, x, d, alpha, nfree) ;
            asa_g (gtemp, xtemp, Com) ;
            dphi = asa_dot (gtemp, d, nfree) ;
            if ( dphi == dphi ) break ;
        }
        if ( i == Parm->nexpand )
        {
            status = -2 ;
            goto Exit ;
        }
        Com->QuadOK = FALSE ;
        rho = Parm->nan_rho ;
    }
    dpsi = dphi - Com->wolfe_hi ;

    /*Find initial interval [a,b] such that dphia < 0, dphib >= 0,
         and phia <= phi0 + feps*Ck */

    a = ZERO ;
    dpsia = dphi0 - Com->wolfe_hi ;
    ngrow = 0 ;
    nshrink = 0 ;
    while ( dpsi < ZERO )
    {
        phi = asa_f (xtemp, Com) ;
        psi = phi - alpha*Com->wolfe_hi ;

        /* if quadstep in effect and quadratic conditions hold,
           check Wolfe condition*/

        if ( Com->QuadOK )
        {
            if ( ngrow == 0 ) fquad = MIN (phi, Com->f0) ;
            if ( phi <= fquad )
            {
                if ( Parm->PrintLevel >= 4 )
                {
                    printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n",
                            alpha, phi, fquad) ;
                }
                if ( asa_Wolfe (alpha, phi, dphi, Com) )
                {
                    status = 0 ;
                    goto Exit ;
                }
            }
        }
        if ( psi <= Com->fpert )
        {
            a = alpha ;
            dpsia = dphi ;
        }
        else
        {
            /* contraction phase, only break at termination or Secant step */
            b = alpha ;
            while ( TRUE )
            {
                alpha = .5*(a+b) ;
                nshrink++ ;
                if ( nshrink > Parm->nexpand )
                {
                    status = 6 ;
                    goto Exit ;
                }
                asa_step (xtemp, x, d, alpha, nfree) ;
                asa_g (gtemp, xtemp, Com) ;
                dphi = asa_dot (gtemp, d, nfree) ;
                dpsi = dphi - Com->wolfe_hi ;
                if ( dpsi >= ZERO ) goto Secant ;
                phi = asa_f (xtemp, Com) ;
                psi = phi - alpha*Com->wolfe_hi ;
                if ( Parm->PrintLevel >= 4 )
                {
                    printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: "
                            "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ;
                }
                if ( Com->QuadOK && (phi <= fquad) )
                {
                    if ( asa_Wolfe (alpha, phi, dphi, Com) )
                    {
                        status = 0 ;
                        goto Exit ;
                    }
                }
                if ( psi <= Com->fpert )
                {
                    a = alpha ;
                    dpsia = dpsi ;
                }
                else
                {
                    b = alpha ;
                }
            }
        }

/* expansion phase */

        ngrow++ ;
        if ( ngrow > Parm->nexpand )
        {
            status = 3 ;
            goto Exit ;
        }
        alphaold = alpha ;
        alpha = MIN (rho*alpha, minstep) ;
        if ( alpha != alphaold )
        {
            asa_step (xtemp, x, d, alpha, nfree) ;
            asa_g (gtemp, xtemp, Com) ;
            dphi = asa_dot (gtemp, d, nfree) ;
            dpsi = dphi - Com->wolfe_hi ;
            if ( Parm->PrintLevel >= 4 )
            {
                printf ("expand,   a: %14.6e alpha: %14.6e phi: "
                         "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ;
            }
        }
        else /* a new constraint is active */
        {
            do /* while statement */
            {
                alphaold = alpha ;
                phiold = phi ;
                if ( alpha < maxstep )
                {
                    alpha = rho*alphaold ;
                    asa_project (xtemp, x, d, alpha, Com) ;
                    phi = asa_f (xtemp, Com) ;
                }
            } while ( phi < phiold ) ;
            if ( alphaold == minstep )
            {
                asa_step (xtemp, x, d, minstep, nfree) ;
                status = -3 ;
            }
            else
            {
                asa_project (xtemp, x, d, alphaold, Com) ;
                asa_g (gtemp, xtemp, Com) ;
                status = -4 ;
            }
            phi = phiold ;
            goto Exit ;