asa_cg.cpp

数学计算程序

            beta = (ykgk - TWO*dphi*ykyk/dkyk)/dkyk ;
/*
    faster: initialize dnorm2 = gnorm2 at start, then
            dnorm2 = gnorm2 + beta**2*dnorm2 - 2.*beta*dphi
            gnorm2 = ||g_{k+1}||^2
            dnorm2 = ||d_{k+1}||^2
            dpi = g_{k+1}' d_k */

            t = -ONE/sqrt (dnorm2*MIN (eta_sq, gnorm2)) ;
            beta = MAX (beta, t) ;

/*    update search direction d = -g + beta*dold */

            gnorm2 = ZERO ;
            for (i = 0; i < n5; i++)
            {
                t = g [i] ;
                d [i] = -t + beta*d [i] ;
                gnorm2 += t*t ;
            }
            for (; i < nfree; )
            {
                t1 = g [i] ;
                d [i] = -t1 + beta*d [i] ;
                i++ ;

                t2 = g [i] ;
                d [i] = -t2 + beta*d [i] ;
                i++ ;

                t3 = g [i] ;
                d [i] = -t3 + beta*d [i] ;
                i++ ;

                t4 = g [i] ;
                d [i] = -t4 + beta*d [i] ;
                i++ ;

                t5 = g [i] ;
                d [i] = -t5 + beta*d [i] ;
                i++ ;

                gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ;
            }
            dphi0 = -gnorm2 + beta*dphi ;
            if ( Parm->debug ) /* Check the dphi0 = d'g */
            {
                t = ZERO ;
                for (j=0; j<nfree; j++)  t = t + d[j]*g[j] ;
                if ( fabs(t-dphi0) > Parm->debugtol*fabs(dphi0) )
                {
                    printf("Warning, dphi0 != d'g!\n");
                    printf("dphi0:%14.6e, d'g:%14.6e\n",dphi0, t) ;
                }
            }
        }
        else
        {
            /* search direction d = -g */
            if ( Parm->PrintLevel >= 3 ) printf ("RESTART CG\n") ;
            ginorm = ZERO ;
            pgnorm = ZERO ;
            gnorm2 = ZERO ;
            for (j = 0; j < nfree; j++)
            {
                xj = xtemp [j] ;
                gj = gtemp [j] ;
                d [j] = -gj ;
                ginorm = MAX (ginorm, fabs (gj)) ;
                gnorm2 += gj*gj ;
                xg = xj - gj ;
                if      ( xg > hi [j] ) xp = hi [j] - xj ;
                else if ( xg < lo [j] ) xp = xj - lo [j] ;
                else                    xp = fabs (gj) ;
                pgnorm = MAX (pgnorm, xp) ;
            }
            for (; j < n; j++)
            {
                xj = x [j] ;
                gj = gtemp [j] ;
                xg = xj - gj ;
                if      ( xg > hi [j] ) xp = hi [j] - xj ;
                else if ( xg < lo [j] ) xp = xj - lo [j] ;
                else                    xp = fabs (gj) ;
                pgnorm = MAX (pgnorm, xp) ;
            }
            if ( asa_tol (pgnorm, Com) )
            {
                status = 0 ;
                for (j = 0; j < nfree; j++)
                {
                    xj = xtemp [j] ;
                    x [j] = xj ;
                    gj = gtemp [j] ;
                    xg = xj - gj ;
                    g [j] = gj ;
                    if      ( xg > hi [j] ) pg [j] = hi [j] - xj ;
                    else if ( xg < lo [j] ) pg [j] = lo [j] - xj ;
                    else                    pg [j] = -gj ;
                }
                for (; j < n; j++)
                {
                    xj = x [j] ;
                    gj = gtemp [j] ;
                    xg = xj - gj ;
                    g [j] = gj ;
                    if      ( xg > hi [j] ) pg [j] = hi [j] - xj ;
                    else if ( xg < lo [j] ) pg [j] = lo [j] - xj ;
                    else                    pg [j] = -gj ;
                }
                goto Exit1 ;
            }
            if ( ginorm < pgnorm*Com->tau2 ) status = -5 ;
            if ( status < -2 )
            {
                sts = ZERO ;
                sty = ZERO ;
                for (j = 0; j < nfree; j++)
                {
                    t = xtemp[j] ;
                    sk = t - x [j] ;
                    x [j] = t ;
                    sts += sk*sk ;

                    t = gtemp [j] ;
                    sty += sk*(t-g [j]) ;
                    g [j] = t ;
                }
                Com->sts = sts ;
                Com->sty = sty ;
                goto Exit ;
            }

            dphi0 = -gnorm2 ;
            asa_copy (x, xtemp, nfree) ;
            asa_copy (g, gtemp, nfree) ;
        }
        if ( !Com->AWolfe )
        {
            if ( fabs (f-Com->f0) <= Parm->AWolfeFac*Ck ) Com->AWolfe = TRUE ;
        }

        if ( Parm->PrintLevel >= 2 )
        {
            printf ("iter: %5i f = %14.6e pgnorm = %14.6e ginorm = %14.6e\n\n",
              (int) iter, f, pgnorm, ginorm) ;
        }

        if ( Parm->debug )
        {
            if ( f > Com->f0 + Ck*Parm->debugtol )
            {
                status = 9 ;
                goto Exit ;
            }
        }

        if ( dphi0 > ZERO )
        {
           status = 5 ;
           goto Exit ;
        }
    }
    status = 2 ;

Exit:
    if ( status < -2 )
    {
        for (j = nfree; j < n; j++) g [j] = gtemp [j] ;
    }
    else
    {
        pgnorm = ZERO ;
        for (j = 0; j < n; j++)
        {
            xj = xtemp [j] ;
            x [j] = xj ;
            gj = gtemp [j] ;
            g [j] = gj ;
            xg = xj - gj ;
            if      ( xg > hi [j] ) xp = hi [j] - xj ;
            else if ( xg < lo [j] ) xp = lo [j] - xj ;
            else                    xp = -gj ;
            pgnorm = MAX (pgnorm, fabs (xp)) ;
            pg [j] = xp ;
        }
    }

Exit1:
    Com->pgnorm = pgnorm ;
    Com->ginorm = ginorm ;
    Com->f = f ;
    Com->f_debug = f ;
    Com->cgfunc += Com->nf - nf ;
    Com->cggrad += Com->ng - ng ;
    Com->cgiter += iter ;
    if ( Parm->PrintLevel >= 2 )
    {
        printf ("iter: %5i f = %14.6e pgnorm = %14.6e ginorm = %14.6e\n\n",
                (int) iter, f, pgnorm, ginorm) ;
    }
    if ( Parm->PrintLevel >= 1 )
    {
        printf ("\nCG Termination status: %i\n", status) ;
        if ( status == -5 )
        {
            printf ("ginorm < tau2*pgnorm without hitting boundary\n") ;
        }
        if ( status == -4 )
        {
            printf ("ginorm >= tau2*pgnorm, many x components hit boundary\n") ;
        }
        else if ( status == -3 )
        {
            printf ("ginorm >= tau2*pgnorm, one x component hits boundary\n") ;
        }
        printf ("proj gradient max norm: %13.6e\n", pgnorm) ;
        printf ("function value:         %13.6e\n", f) ;
        printf ("cg iterations:          %13.6e\n", (double) iter) ;
        printf ("function evaluations:   %13.6e\n", (double) Com->nf - nf) ;
        printf ("gradient evaluations:   %13.6e\n", (double) Com->ng - ng) ;
    }
    return (status) ;
}

/* =========================================================================
   === asa_Wolfe ===========================================================
   =========================================================================
   Check whether the Wolfe or the approximate Wolfe conditions are satisfied
   ========================================================================= */
int asa_Wolfe
(
    double       alpha , /* stepsize */
    double           f , /* function value associated with stepsize alpha */
    double        dphi , /* derivative value associated with stepsize alpha */
    asa_com        *Com  /* cg com */
)
{
    if ( dphi >= Com->wolfe_lo )
    {

        /* test original Wolfe conditions */
        if ( f - Com->f0 <= alpha*Com->wolfe_hi )
        {
            if ( Com->cgParm->PrintLevel >= 4 )
            {
                printf ("wolfe f: %14.6e f0: %14.6e dphi: %14.6e\n",
                         f, Com->f0, dphi) ;
            }
            return (1) ;
        }
        /* test approximate Wolfe conditions */
        else if ( Com->AWolfe )
        {
            if ( (f <= Com->fpert) && (dphi <= Com->awolfe_hi) )
            {
                if ( Com->cgParm->PrintLevel >= 4 )
                {
                    printf ("f: %14.6e fpert: %14.6e dphi: %14.6e awolf_hi: "
                            "%14.6e\n", f, Com->fpert, dphi, Com->awolfe_hi) ;
                }
                return (1) ;
            }
        }
    }
    return (0) ;
}

/* =========================================================================
   === asa_tol =============================================================
   =========================================================================
   Check for convergence
   ========================================================================= */
int asa_tol
(
    double      pgnorm, /* projected gradient sup-norm */
    asa_com       *Com
)
{
    /*StopRule = T => |grad|_infty <=max (tol, |grad|_infty*StopFact)
                 F => |grad|_infty <= tol*(1+|f|)) */
    if ( Com->asaParm->StopRule )
    {
        if ( pgnorm <= Com->tol ) return (1) ;
    }
    else if ( pgnorm <= Com->tol*(ONE + fabs (Com->f)) ) return (1) ;
    return (0) ;
}

/* =========================================================================
   === asa_step ============================================================
   =========================================================================
   Compute xtemp = x + alpha d
   ========================================================================= */
void asa_step
(
    double *xtemp , /*output vector */
    double     *x , /* initial vector */
    double     *d , /* search direction */
    double  alpha , /* stepsize */
    INT32         n   /* length of the vectors */
)
{
    INT32 n5, i ;
    n5 = n % 5 ;
    for (i = 0; i < n5; i++) xtemp [i] = x[i] + alpha*d[i] ;
    for (; i < n; i += 5)
    {
        xtemp [i]   = x [i]   + alpha*d [i] ;
        xtemp [i+1] = x [i+1] + alpha*d [i+1] ;
        xtemp [i+2] = x [i+2] + alpha*d [i+2] ;
        xtemp [i+3] = x [i+3] + alpha*d [i+3] ;
        xtemp [i+4] = x [i+4] + alpha*d [i+4] ;
    }
}

/* =========================================================================
   === asa_line ============================================================
   =========================================================================
   Approximate Wolfe line search routine
   ========================================================================= */
int asa_line
(
    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, dphia, b, dphib, c, alpha, phi, dphi, alphaold, phiold,
           a0, da0, b0, db0, width, fquad, rho, minstep, maxstep,
           *x, *xtemp, *d, *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 ;
    }

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

    a = ZERO ;
    dphia = dphi0  ;
    ngrow = 0 ;
    nshrink = 0 ;
    while ( dphi < ZERO )
    {
        phi = asa_f (xtemp, Com) ;

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

        if ( Com->QuadOK )
        {
            if ( ngrow == 0 ) fquad = MIN (phi, Com->f0) ;