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) ;