asa_cg.cpp

数学计算程序

   ========================================================================= */
int asa_descent /*  return:
                      -5 (ginorm < tau2*pgnorm without hitting boundary)
                      -4 (ginorm >=tau2*pgnorm, many x components hit boundary)
                      -3 (ginorm >=tau2*pgnorm, one x component hits boundary)
                      -2 (function value became nan)
                      -1 (starting function value is nan)
                       0 (convergence tolerance satisfied)
                       1 (change in func <= feps*|f|)
                       2 (total iterations exceeded maxit)
                       3 (slope always negative in line search)
                       4 (number secant iterations exceed nsecant)
                       5 (search direction not a descent direction)
                       6 (line search fails in initial interval)
                       7 (line search fails during bisection)
                       8 (line search fails during interval update)
                       9 (debugger is on and the function value increases)*/
(
    asa_com *Com
)
{
    int     i, iter, j, maxit, n, n5, nfree, nf, ng, nrestart,
            status ;
    double  delta2, eta_sq, Qk, Ck, pgnorm, ginorm,
            f, ftemp, gnorm, xnorm, gnorm2, dnorm2, denom,
            t, t1, t2, t3, t4, t5, dphi, dphi0, alpha, talpha,
            xj, gj, xg, xp, sts, sty, sk,
            yk, ykyk, ykgk, dkyk, yk1, yk2, yk3, yk4, yk5, beta,
            *x, *d, *g, *xtemp, *gtemp, *lo, *hi, *pg ;

    asacg_parm *Parm ;
    asa_parm *asaParm ;

/* initialization */

    x = Com->x ;
    lo = Com->lo ;
    hi = Com->hi ;
    n = Com->n ;
    d = Com->d ;
    g = Com->g ;
    xtemp = Com->xtemp ;
    gtemp = Com->gtemp ;
    pg = Com->pg ;
    nfree = Com->nfree ;
    nf = Com->nf ;
    ng = Com->ng ;
    pgnorm = Com->pgnorm ;
    ginorm = Com->ginorm ;
    Parm = Com->cgParm ;
    asaParm = Com->asaParm ;

    if ( Parm->PrintLevel >= 1 )
    {
        printf ("Dimension in CG, nfree = %i\n", nfree) ;
    }

    /* the conjugate gradient algorithm is restarted every nrestart iteration */
    nrestart = (INT32) (((double) nfree)*Parm->restart_fac) ;

    /* abort when number of iterations reaches maxit in one pass through cg */
    if ( Parm->maxit_fac == INF ) Com->cgmaxit = maxit = INT_INF ;
    else Com->cgmaxit = maxit = (INT32) (((double) n)*Parm->maxit_fac) ;

    n5 = nfree % 5 ;
    f = Com->f ;

    Ck = ZERO ;
    Qk = ZERO ;

/* initial function and gradient evaluations, initial direction */

    Com->f0 = f + f ;
    xnorm = asa_max (x, nfree) ;
    gnorm = ZERO ;
    gnorm2 = ZERO ;
    for (i = 0; i < n5; i++)
    {
        t = g [i] ;
        d [i] = -t ;
        gnorm2 += t*t ;
        if ( gnorm < fabs (t) ) gnorm = fabs (t) ;
    }
    for (; i < nfree;)
    {
        t1 = g [i] ;
        d [i] = -t1 ;
        if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ;
        i++ ;

        t2 = g [i] ;
        d [i] = -t2 ;
        if ( gnorm < fabs (t2) ) gnorm = fabs (t2) ;
        i++ ;

        t3 = g [i] ;
        d [i] = -t3 ;
        if ( gnorm < fabs (t3) ) gnorm = fabs (t3) ;
        i++ ;

        t4 = g [i] ;
        d [i] = -t4 ;
        if ( gnorm < fabs (t4) ) gnorm = fabs (t4) ;
        i++ ;

        t5 = g [i] ;
        d [i] = -t5 ;
        if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ;
        i++ ;

        gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ;
    }
    /* check that starting function value is nan */
    if ( f != f )
    {
        status = -1 ;
        goto Exit ;
    }

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

    dphi0 = -gnorm2 ;
    delta2 = 2*Parm->delta - ONE ;
    eta_sq = Parm->eta*Parm->eta ;
    alpha = Parm->step ;
    if ( alpha == ZERO )
    {
        alpha = Parm->psi0*xnorm/gnorm ;
        if ( xnorm == ZERO )
        {
            if ( f != ZERO ) alpha = Parm->psi0*fabs (f)/gnorm2 ;
            else             alpha = ONE ;
        }
    }

/*  start the conjugate gradient iteration
    alpha starts as old step, ends as final step for current iteration
    f is function value for alpha = 0
    Com->QuadOK = TRUE means that a quadratic step was taken */

    for (iter = 1; iter <= maxit; iter++)
    {
        Com->QuadOK = FALSE ;
        alpha = Parm->psi2*alpha ;
        asa_maxstep (x, d, Com) ;
        if ( Parm->QuadStep )
        {
            if ( f != ZERO ) t = fabs ((f-Com->f0)/f) ;
            else             t = ONE ;
            if ( t > Parm->QuadCutOff )       /* take provisional step talpha */
            {
                talpha = Parm->psi1*alpha ;
                if ( Com->minstep >= asaParm->parm1*talpha )
                {
                    talpha = MIN (talpha, asaParm->parm2*Com->minstep) ;
                    asa_step (xtemp, x, d, talpha, nfree) ;
                    /*provisional function value*/
                    ftemp = asa_f (xtemp, Com) ;

                    /* check if function value is nan */
                    if ( ftemp != ftemp ) /* reduce stepsize */
                    {
                        for (i = 0; i < Parm->nexpand; i++)
                        {
                            talpha /= Parm->rho ;
                            asa_step (xtemp, x, d, talpha, nfree) ;
                            ftemp = asa_f (xtemp, Com) ;
                            if ( ftemp == ftemp ) break ;
                        }
                        if ( i == Parm->nexpand )
                        {
                            status = -2 ;
                            goto Exit ;
                        }
                    }

                    if ( ftemp < f )              /* check if quadstep > 0 */
                    {
                       denom = TWO*(((ftemp-f)/talpha)-dphi0) ;
                       if ( denom > ZERO )    /* try a quadratic fit step */
                       {
                           Com->QuadOK = TRUE ;
                           alpha = -dphi0*talpha/denom ;
                       }
                    }
                }
            }
        }
        Com->f0 = f ;                          /* f0 saved as prior value */
        if ( Parm->PrintLevel >= 3 )
        {
            printf ("minstep =%14.6e, maxstep =%14.6e \n",
                     Com->minstep, Com->maxstep) ;
            printf ("QuadOK: %2i initial a: %14.6e f0: %14.6e dphi0: %14.6e\n",
                    Com->QuadOK, alpha, Com->f0, dphi0) ;
            if ( (alpha > Com->minstep) && Com->QuadOK )
            {
                printf("Quadratic step > minstep to boundary\n") ;
            }
        }

        /* parameters in Wolfe, approximate Wolfe conditions, and in update */
        Qk = Parm->Qdecay*Qk + ONE ;
        Ck = Ck + (fabs (f) - Ck)/Qk ;        /* average cost magnitude */

        if ( Parm->PertRule ) Com->fpert = f + Parm->eps*Ck ;
        else                  Com->fpert = f + Parm->eps ;

        Com->wolfe_hi = Parm->delta*dphi0 ;
        Com->wolfe_lo = Parm->sigma*dphi0 ;
        Com->awolfe_hi = delta2*dphi0 ;
        Com->alpha = alpha ;/* either double prior step or quadratic fit step */
        Com->f = f ;

        if ( Com->AWolfe )                  /* approximate Wolfe line search*/
        {
            if ( Parm->PrintLevel >= 3 )
            {
                printf ("Perform approximate Wolfe line search\n") ;
            }

            status = asa_line (dphi0, Com) ;
        }
        else                                  /* ordinary Wolfe line search */
        {
            if ( Parm->PrintLevel >= 3 )
            {
                 printf ("Perform ordinary Wolfe line search\n") ;
            }
            status = asa_lineW (dphi0, Com) ;
        }
        /* if ordinary Wolfe line search fails, possibly try approximate
           Wolfe line search*/
        if ( (status > 0) && !Com->AWolfe && (Parm->AWolfeFac > ZERO) )
        {
            Com->AWolfe = TRUE ;
            if ( Parm->PrintLevel >= 3 )
            {
                printf ("Ordinary Wolfe line search fails, "
                        "try approximate Wolfe line search\n") ;
            }

            status = asa_line (dphi0, Com) ;
        }

        alpha = Com->alpha ;
        f = Com->f ;
        dphi = Com->df ;

        if ( (status > 0) || (status == -1) || (status == -2) ) goto Exit ;

        /*Test for convergence to within machine epsilon
          [set feps to zero to remove this test] */

        if ( (-alpha*dphi0 <= Parm->feps*fabs (f)) && (status == 0) )
        {
            status = 1 ;
            goto Exit ;
        }

        /* compute beta, yk2, gnorm, gnorm2, dnorm2, update x and g */
        if ( iter % nrestart != 0 )
        {
            ginorm = ZERO ;
            pgnorm = ZERO ;
            for (j = 0; j < nfree; j++)
            {
                xj = xtemp [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) ;
                ginorm = MAX (ginorm, fabs (gj)) ;
            }
            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 ;
            }

            asa_copy (x, xtemp, nfree) ;
            dnorm2 = ZERO ;
            for (j = 0; j < n5; j++) dnorm2 = dnorm2 + d [j]*d [j] ;
            for (; j < nfree; j += 5)
            {
                dnorm2 = dnorm2 + d [j]*d [j] + d [j+1]*d [j+1]
                                              + d [j+2]*d [j+2]
                                              + d [j+3]*d [j+3]
                                              + d [j+4]*d [j+4] ;
            }
            ykyk = ZERO ;
            ykgk = ZERO ;
            for (j = 0; j < n5; j++)
            {
                t = gtemp [j] ;
                yk = t - g [j] ;
                ykyk += yk*yk ;
                ykgk += yk*t ;
                g [j] = t ;
            }
            for (j = n5; j < nfree; )
            {
                t1 = gtemp [j] ;
                yk1 = t1 - g [j] ;
                g [j] = t1 ;
                j++ ;

                t2 = gtemp [j] ;
                yk2 = t2 - g [j] ;
                g [j] = t2 ;
                j++ ;

                t3 = gtemp [j] ;
                yk3 = t3 - g [j] ;
                g [j] = t3 ;
                j++ ;

                t4 = gtemp [j] ;
                yk4 = t4 - g [j] ;
                g [j] = t4 ;
                j++ ;

                t5 = gtemp [j] ;
                yk5 = t5 - g [j] ;
                g [j] = t5 ;
                j++ ;

                ykyk += yk1*yk1 + yk2*yk2 + yk3*yk3 + yk4*yk4 + yk5*yk5 ;
                ykgk += yk1*t1 + yk2*t2 + yk3*t3 + yk4*t4 + yk5*t5 ;
            }

            dkyk = dphi - dphi0 ;