/* =========================================================================
   ============================ CG_DESCENT =================================
   =========================================================================
       ________________________________________________________________
      |      A conjugate gradient method with guaranteed descent       |
      |             C-code Version 1.1  (October 6, 2005)              |
      |                    Version 1.2  (November 14, 2005)            |
      |                    Version 2.0  (September 23, 2007)           |
      |                    Version 3.0  (May 18, 2008)                 |
      |           William W. Hager    and   Hongchao Zhang             |
      |          hager@math.ufl.edu       hzhang@math.ufl.edu          |
      |                   Department of Mathematics                    |
      |                     University of Florida                      |
      |                 Gainesville, Florida 32611 USA                 |
      |                      352-392-0281 x 244                        |
      |                                                                |
      |                 Copyright by William W. Hager                  |
      |                                                                |
      |          http://www.math.ufl.edu/~hager/papers/CG              |
      |________________________________________________________________|
       ________________________________________________________________
      |This program is free software; you can redistribute it and/or   |
      |modify it under the terms of the GNU General Public License as  |
      |published by the Free Software Foundation; either version 2 of  |
      |the License, or (at your option) any later version.             |
      |This program is distributed in the hope that it will be useful, |
      |but WITHOUT ANY WARRANTY; without even the implied warranty of  |
      |MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the   |
      |GNU General Public License for more details.                    |
      |                                                                |
      |You should have received a copy of the GNU General Public       |
      |License along with this program; if not, write to the Free      |
      |Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
      |MA  02110-1301  USA                                             |
      |________________________________________________________________|*/

#include "cg_user.h"
#include "cg_descent.h"
int cg_descent /*  return:
                      -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)
                      10 (out of memory) */
(
    double            *x, /* input: starting guess, output: the solution */
    INT                n, /* problem dimension */
    cg_stats       *Stat, /* structure with statistics (can be NULL) */
    cg_parameter  *UParm, /* user parameters, NULL = use default parameters */
    double      grad_tol, /* StopRule = 1: |g|_infty <= max (grad_tol,
                                           StopFac*initial |g|_infty) [default]
                             StopRule = 0: |g|_infty <= grad_tol(1+|f|) */
    double      (*value) (double *, INT),  /* f = value (x, n) */
    void         (*grad) (double *, double *, INT), /* grad (g, x, n) */
    double    (*valgrad) (double *, double *, INT), /* f = valgrad (g, x, n),
                          NULL = compute value & gradient using value & grad */
    double         *Work  /* either size 4n work array or NULL */
)
{
    INT     n5, iter, maxit, nrestart, i ;
    int     status, StopRule ;
    double  delta2, eta_sq, Qk, Ck,
            f, ftemp, gnorm, xnorm, gnorm2, dnorm2, denom,
            t, t1, t2, t3, t4, t5, dphi, dphi0, alpha, talpha,
            yk, ykyk, ykgk, dkyk, yk1, yk2, yk3, yk4, yk5, beta, tol,
           *d, *g, *xtemp, *gtemp, *work ;
    cg_parameter *Parm, ParmStruc ;
    cg_com Com ;

/* initialize the parameters */

    if ( UParm == NULL )
    {
        Parm = &ParmStruc ;
        cg_default (Parm) ;
    }
    else Parm = UParm ;
    Com.Parm = Parm ;

    if ( Parm->PrintParms ) cg_printParms (Parm) ;

    /* allocate work arrays */
    if ( Work == NULL ) work = malloc (4*n*sizeof (double)) ;
    else                work = Work ;
    if ( work == NULL )
    {
        printf ("Insufficient memory for specified problem dimension %e\n",
                 (double) n) ;
        status = 10 ;
        return (status) ;
    }
    Com.x = x ;
    Com.d = d = work ;
            g = d+n ;
    Com.xtemp = xtemp = g+n ;
    Com.gtemp = gtemp = xtemp+n ;
    Com.n = n ;          /* problem dimension */
    Com.nf = (INT) 0 ;   /* number of function evaluations */
    Com.ng = (INT) 0 ;   /* number of gradient evaluations */
    Com.AWolfe = Parm->AWolfe ; /* do not touch user's AWolfe */
    Com.cg_value = value ;
    Com.cg_grad = grad ;
    Com.cg_valgrad = valgrad ;
    StopRule = Parm->StopRule ;

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

    /* abort when number of iterations reaches maxit */
    if ( Parm->maxit_fac == INF ) maxit = INT_INF ;
    else                          maxit = (INT) (((double) n)*Parm->maxit_fac) ;
    
    f = ZERO ;
    n5 = n % 5 ;

    Ck = ZERO ;
    Qk = ZERO ;

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

    f = cg_fg (g, x, &Com) ;
    Com.f0 = f + f ;
    xnorm = ZERO ;
    for (i = 0; i < n5; i++) if ( xnorm < fabs (x [i]) ) xnorm = fabs (x [i]) ;
    for (; i < n; i += 5)
    {
         if ( xnorm < fabs (x [i]  ) ) xnorm = fabs (x [i]  ) ;
         if ( xnorm < fabs (x [i+1]) ) xnorm = fabs (x [i+1]) ;
         if ( xnorm < fabs (x [i+2]) ) xnorm = fabs (x [i+2]) ;
         if ( xnorm < fabs (x [i+3]) ) xnorm = fabs (x [i+3]) ;
         if ( xnorm < fabs (x [i+4]) ) xnorm = fabs (x [i+4]) ;
    }
    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 < n;)
    {
        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->StopRule ) tol = MAX (gnorm*Parm->StopFac, grad_tol) ;
    else                  tol = grad_tol ;

    if ( Parm->PrintLevel >= 1 )
    {
        printf ("iter: %5i f = %14.6e gnorm = %14.6e AWolfe = %2i\n",
          (int) 0, f, gnorm, Com.AWolfe) ;
    }

    if ( cg_tol (f, gnorm, StopRule, tol) )
    {
        status = 0 ;
        goto Exit ;
    }

    dphi0 = -gnorm2 ;
    delta2 = 2*Parm->delta - ONE ;
    eta_sq = Parm->eta*Parm->eta ;
    alpha = Parm->step ;
    if ( alpha == 0. )
    {
        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
    QuadOK = TRUE means that a quadratic step was taken */
 
    for (iter = 1; iter <= maxit; iter++)
    {
        Com.QuadOK = FALSE ;
        alpha = Parm->psi2*alpha ;
        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 ;
                cg_step (xtemp, x, d, talpha, n) ;
                ftemp = cg_f (xtemp, &Com) ;  /* provisional function value */

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

                if ( ftemp < f )              /* check if quadstep > 0 */
                {
                   denom = 2.*(((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 >= 1 )
        {
            printf ("QuadOK: %2i initial a: %14.6e f0: %14.6e dphi: %14.6e\n",
                    Com.QuadOK, alpha, Com.f0, dphi0) ;
        }

/* parameters in Wolfe and 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 prior step or quadratic fit step */
        Com.f = f ;
        
        if ( Com.AWolfe ) status = cg_line (dphi0, &Com) ; /* approx. Wolfe */
        else              status = cg_lineW (dphi0, &Com) ;/* ordinary Wolfe */
        if ( (status > 0) && !Com.AWolfe )/*try approximate Wolfe line search*/
        {
            if ( Parm->PrintLevel >= 1 )
            {
                 printf ("\nWOLFE LINE SEARCH FAILS\n") ;
            }
            Com.AWolfe = TRUE ;
            status = cg_line (dphi0, &Com) ;
        }

        alpha = Com.alpha ;
        f = Com.f ;
        dphi = Com.df ;

        if ( status ) 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 = 1 ;
            goto Exit ;
        }

/* compute beta, yk2, gnorm, gnorm2, dnorm2, update x and g */

        if ( iter % nrestart != 0 )
        {
            cg_copy (x, xtemp, n) ;
            dnorm2 = ZERO ;
            for (i = 0; i < n5; i++) dnorm2 = dnorm2 + d [i]*d [i] ;
            for (; i < n; i += 5)
            {
                dnorm2 = dnorm2 + d [i]*d [i] + d [i+1]*d [i+1]
                                              + d [i+2]*d [i+2]
                                              + d [i+3]*d [i+3]
                                              + d [i+4]*d [i+4] ;
            }
            gnorm = ZERO ;
            ykyk = ZERO ;
            ykgk = ZERO ;
            for (i = 0; i < n5; i++)
            {
                t = gtemp [i] ;
                if ( gnorm < fabs (t) ) gnorm = fabs (t) ;
                yk = t - g [i] ;
                g [i] = t ;
                ykgk += yk*t ;
                ykyk += yk*yk ;
            }
            for (; i < n; )
            {
                t1 = gtemp [i] ;
                yk1 = t1 - g [i] ;
                g [i] = t1 ;
                if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ;
                i++ ;

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

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

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

                t5 = gtemp [i] ;
                yk5 = t5 - g [i] ;
                g [i] = t5 ;
                if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ;

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

            if ( cg_tol (f, gnorm, StopRule, tol) )
            {
                status = 0 ;
                goto Exit ;
            }
            dkyk = dphi - dphi0 ;
            beta = (ykgk - 2.*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 < n; )
            {
                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 (i = 0; i < n; i++)  t = t + d [i]*g [i] ;
                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 >= 1 ) printf ("RESTART CG\n") ;
            gnorm = ZERO ;
            gnorm2 = ZERO ;
            cg_copy (x, xtemp, n) ;
            for (i = 0; i < n5; i++)
            {
                t = gtemp [i] ;
                g [i] = t ;
                d [i] = -t ;
                if ( gnorm < fabs (t) ) gnorm = fabs (t) ;
                gnorm2 += t*t ;
            }
            for (; i < n; )
            {
                t1 = gtemp [i] ;
                g [i] = t1 ;
                d [i] = -t1 ;
                if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ;
                i++ ;

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

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

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

                t5 = gtemp [i] ;
                g [i] = t5 ;
                d [i] = -t5 ;
                if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ;
                i++ ;
                gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ;
            }
            if ( cg_tol (f, gnorm, StopRule, tol) )
            {
                status = 0 ;
                goto Exit ;
            }
            dphi0 = -gnorm2 ;
        }
        if ( !Com.AWolfe )
        {
            if ( fabs (f-Com.f0) < Parm->AWolfeFac*Ck ) Com.AWolfe = TRUE ;
        }
    
        if ( Parm->PrintLevel >= 1 )
        {
            printf ("\niter: %5i f = %14.6e gnorm = %14.6e AWolfe = %2i\n",
               (int) iter, f, gnorm, Com.AWolfe) ;
        }

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

Exit:
    if ( Stat != NULL )
    {
        Stat->f = f ;
        Stat->gnorm = gnorm ;
        Stat->nfunc = Com.nf ;
        Stat->ngrad = Com.ng ;
        Stat->iter = iter ;
    }
    if ( status > 2 )
    {
        gnorm = ZERO ;
        for (i = 0; i < n; i++)
        {
            x [i] = xtemp [i] ;
            g [i] = gtemp [i] ;
            t = fabs (g [i]) ;
            gnorm = MAX (gnorm, t) ;
        }
        if ( Stat != NULL ) Stat->gnorm = gnorm ;
    }
    if ( Parm->PrintFinal || Parm->PrintLevel >= 1 )
    {
        const char mess1 [] = "Possible causes of this error message:" ;
        const char mess2 [] = "   - your tolerance may be too strict: "
                              "grad_tol = " ;
        const char mess3 [] = "Line search fails" ;
        const char mess4 [] = "   - your gradient routine has an error" ;
        const char mess5 [] = "   - the parameter epsilon in cg_descent_c.parm "
                              "is too small" ;
        printf ("\nTermination status: %i\n", status) ;
        if ( status == -2 )
        {
            printf ("At iteration %10.0f function value became nan\n",
                    (double) iter) ;
        }
        else if ( status == -1 )
        {
            printf ("Objective function value is nan at starting point\n") ;
        }
        else if ( status == 0 )
        {
            printf ("Convergence tolerance for gradient satisfied\n") ;
        }
        else if ( status == 1 )
        {
            printf ("Terminating since change in function value "
                    "<= feps*|f|\n") ;
        }
        else if ( status == 2 )
        {
            printf ("Number of iterations exceed specified limit\n") ;
            printf ("Iterations: %10.0f maxit: %10.0f\n",
                    (double) iter, (double) maxit) ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, grad_tol) ;
        }
        else if ( status == 3 )
        {
            printf ("Slope always negative in line search\n") ;
            printf ("%s\n", mess1) ;
            printf ("   - your cost function has an error\n") ;
            printf ("%s\n", mess4) ;
        }
        else if ( status == 4 )
        {
            printf ("Line search fails, too many secant steps\n") ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, grad_tol) ;
        }
        else if ( status == 5 )
        {
            printf ("Search direction not a descent direction\n") ;
        }
        else if ( status == 6 ) /* line search fails */
        {
            printf ("%s\n", mess3) ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, grad_tol) ;
            printf ("%s\n", mess4) ;
            printf ("%s\n", mess5) ;
        }
        else if ( status == 7 ) /* line search fails */
        {
            printf ("%s\n", mess3) ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, grad_tol) ;
        }
        else if ( status == 8 ) /* line search fails */
        {
            printf ("%s\n", mess3) ;
            printf ("%s\n", mess1) ;
            printf ("%s %e\n", mess2, grad_tol) ;
            printf ("%s\n", mess4) ;
            printf ("%s\n", mess5) ;
        }
        else if ( status == 9 )
        {
            printf ("Debugger is on, function value does not improve\n") ;
            printf ("new value: %25.16e old value: %25.16e\n", f, Com.f0) ;
        }
        else if ( status == 10 )
        {
            printf ("Insufficient memory\n") ;
        }

        printf ("maximum norm for gradient: %13.6e\n", gnorm) ;
        printf ("function value:            %13.6e\n\n", f) ;
        printf ("cg  iterations:          %10.0f\n", (double) iter) ;
        printf ("function evaluations:    %10.0f\n", (double) Com.nf) ;
        printf ("gradient evaluations:    %10.0f\n", (double) Com.ng) ;
        printf ("===================================\n\n") ;
    }
    if ( Work == NULL ) free (work) ;
    return (status) ;
}

/* =========================================================================
   === cg_default ==========================================================
   =========================================================================
   Set default conjugate gradient parameter values. If the parameter argument
   of cg_descent is NULL, this routine is called by cg_descent automatically.
   If the user wishes to set parameter values, then the cg_parameter structure
   should be allocated in the main program. The user could call cg_default
   to initialize the structure, and then individual elements in the structure
   could be changed, before passing the structure to cg_descent.
   =========================================================================*/
void cg_default
(
    cg_parameter   *Parm
)
{
    /* T => print final function value
       F => no printout of final function value */
    Parm->PrintFinal = TRUE ;

   /* Level 0 = no printing, ... , Level 3 = maximum printing */
    Parm->PrintLevel = 0 ;

    /* T => print parameters values
       F => do not display parmeter values */
    Parm->PrintParms = FALSE ;

    /* T => use approximate Wolfe line search
       F => use ordinary Wolfe line search, switch to approximate Wolfe when
                |f_k+1-f_k| < AWolfeFac*C_k, C_k = average size of cost */
    Parm->AWolfe = FALSE ;
    Parm->AWolfeFac = 1.e-3 ;

    /* factor in [0, 1] used to compute average cost magnitude C_k as follows:
       Q_k = 1 + (Qdecay)Q_k-1, Q_0 = 0,  C_k = C_k-1 + (|f_k| - C_k-1)/Q_k */
    Parm->Qdecay = .7 ;

    /* Stop Rules:
       T => ||grad||_infty <= max(grad_tol, initial |grad|_infty*StopFact)
       F => ||grad||_infty <= grad_tol*(1 + |f_k|) */
    Parm->StopRule = TRUE ;
    Parm->StopFac = 0.e-12 ;

    /* T => estimated error in function value is eps*Ck,
       F => estimated error in function value is eps */
    Parm->PertRule = TRUE ;
    Parm->eps = 1.e-6 ;

    /* T => attempt quadratic interpolation in line search when
                |f_k+1 - f_k|/f_k <= QuadCutoff
       F => no quadratic interpolation step */
    Parm->QuadStep = TRUE ;
    Parm->QuadCutOff = 1.e-12 ;

    /* T => check that f_k+1 - f_k <= debugtol*C_k
       F => no checking of function values */
    Parm->debug = FALSE ;
    Parm->debugtol = 1.e-10 ;

    /* if step is nonzero, it is the initial step of the initial line search */
    Parm->step = ZERO ;

    /* abort cg after maxit_fac*n iterations */
    Parm->maxit_fac = INF ;

    /* maximum number of times the bracketing interval grows or shrinks
       in the line search is nexpand */
    Parm->nexpand = (int) 50 ;

    /* maximum number of secant iterations in line search is nsecant */
    Parm->nsecant = (int) 50 ;

    /* conjugate gradient method restarts after (n*restart_fac) iterations */
    Parm->restart_fac = ONE ;

    /* stop when -alpha*dphi0 (estimated change in function value) <= feps*|f|*/
    Parm->feps = ZERO ;

    /* after encountering nan, growth factor when searching for
       a bracketing interval */
    Parm->nan_rho = 1.3 ;

    /* Wolfe line search parameter, range [0, .5]
       phi (a) - phi (0) <= delta phi'(0) */
    Parm->delta = .1 ;

    /* Wolfe line search parameter, range [delta, 1]
       phi' (a) >= sigma phi' (0) */
    Parm->sigma = .9 ;

    /* decay factor for bracket interval width in line search, range (0, 1) */
    Parm->gamma = .66 ;

    /* growth factor in search for initial bracket interval */
    Parm->rho = 5. ;

    /* conjugate gradient parameter beta_k must be >= eta*||d_k||_2 */
    Parm->eta = .01 ;

    /* starting guess for line search =
         psi0 ||x_0||_infty over ||g_0||_infty if x_0 != 0
         psi0 |f(x_0)|/||g_0||_2               otherwise */
    Parm->psi0 = .01 ;      /* factor used in starting guess for iteration 1 */

    /* for a QuadStep, function evalutated at psi1*previous step */
    Parm->psi1 = .1 ;

    /* when starting a new cg iteration, our initial guess for the line
       search stepsize is psi2*previous step */
    Parm->psi2 = 2. ;
}

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

/* test original Wolfe conditions */

        if ( f - Com->f0 <= alpha*Com->wolfe_hi )
        {
            if ( Com->Parm->PrintLevel >= 2 )
            {
                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->Parm->PrintLevel >= 2 )
                {
                    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) ;
}

/* =========================================================================
   ==== cg_f ===============================================================
   Evaluate the function
   =========================================================================*/
double cg_f
(
    double   *x,
    cg_com *Com
)
{
    double f ;
    f = Com->cg_value (x, Com->n) ;
    Com->nf++ ;
    return (f) ;
}

/* =========================================================================
   ==== cg_g ===============================================================
   Evaluate the gradient
   =========================================================================*/
void cg_g
(
    double   *g,
    double   *x,
    cg_com *Com
)
{
    Com->cg_grad (g, x, Com->n) ;
    Com->ng++ ;
}

/* =========================================================================
   ==== cg_fg ==============================================================
   Evaluate the function and gradient
   =========================================================================*/
double cg_fg
(
    double   *g,
    double   *x,
    cg_com *Com
)
{
    double f ;
    if ( Com->cg_valgrad != NULL ) f = Com->cg_valgrad (g, x, Com->n) ;
    else
    {
        Com->cg_grad (g, x, Com->n) ;
        f = Com->cg_value (x, Com->n) ;
    }
    Com->nf++ ;
    Com->ng++ ;
    return (f) ;
}

/* =========================================================================
   ==== cg_tol =============================================================
   =========================================================================
   Check for convergence
   ========================================================================= */
int cg_tol
(
    double         f, /* function value associated with stepsize */
    double     gnorm, /* gradient sup-norm */
    int     StopRule, /* T => |grad|_infty <=max (tol, |grad|_infty*StopFact)
                         F => |grad|_infty <= tol*(1+|f|)) */
    double       tol  /* tolerance */
)
{
    if ( StopRule )
    {
        if ( gnorm <= tol ) return (1) ;
    }
    else if ( gnorm <= tol*(ONE + fabs (f)) ) return (1) ;
    return (0) ;
}

/* =========================================================================
   ==== cg_dot =============================================================
   =========================================================================
   Compute dot product of x and y, vectors of length n
   ========================================================================= */
double cg_dot
(
    double *x, /* first vector */
    double *y, /* second vector */
    INT     n /* length of vectors */
)
{
    INT i, n5 ;
    double t ;
    t = 0. ;
    n5 = n % 5 ;
    for (i = 0; i < n5; i++) t += x [i]*y [i] ;
    for (; i < n; i += 5)
    {
        t += x [i]*y[i] + x [i+1]*y [i+1] + x [i+2]*y [i+2]
                        + x [i+3]*y [i+3] + x [i+4]*y [i+4] ;
    }
    return (t) ;
}

/* =========================================================================
   === cg_copy =============================================================
   =========================================================================
   Copy vector x into vector y
   ========================================================================= */
void cg_copy
(
    double *y, /* output of copy */
    double *x, /* input of copy */
    int     n  /* length of vectors */
)
{
    int j, n10 ;
    n10 = n % 10 ;
    for (j = 0; j < n10; j++) y [j] = x [j] ;
    for (; j < n; j += 10)
    {
        y [j] = x [j] ;
        y [j+1] = x [j+1] ;
        y [j+2] = x [j+2] ;
        y [j+3] = x [j+3] ;
        y [j+4] = x [j+4] ;
        y [j+5] = x [j+5] ;
        y [j+6] = x [j+6] ;
        y [j+7] = x [j+7] ;
        y [j+8] = x [j+8] ;
        y [j+9] = x [j+9] ;
    }
}

/* =========================================================================
   ==== cg_step ============================================================
   =========================================================================
   Compute xtemp = x + alpha d
   ========================================================================= */
void cg_step
(
    double *xtemp, /*output vector */
    double     *x, /* initial vector */
    double     *d, /* search direction */
    double  alpha, /* stepsize */
    INT         n  /* length of the vectors */
)
{
    INT 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] ;
    }
}

/* =========================================================================
   ==== cg_line ============================================================
   =========================================================================
   Approximate Wolfe line search routine
   ========================================================================= */
int cg_line
(
    double  dphi0, /* function derivative at starting point (alpha = 0) */
    cg_com   *Com  /* cg com structure */
)
{
    INT n, iter ;
    int i, nsecant, nshrink, ngrow, status ;
    double a, dphia, b, dphib, c, alpha, phi, dphi,
           a0, da0, b0, db0, width, fquad, rho, *x, *xtemp, *d, *gtemp ;
    cg_parameter *Parm ;

    Parm = Com->Parm ;
    if ( Parm->PrintLevel >= 1 ) printf ("Approximate Wolfe line search\n") ;
    alpha = Com->alpha ;
    phi = Com->f ;
    n = Com->n ;
    x = Com->x ;         /* current iterate */
    xtemp = Com->xtemp ; /* x + alpha*d */
    d = Com->d ;         /* current search direction */
    gtemp = Com->gtemp ; /* gradient at x + alpha*d */
    rho = Parm->rho ;
    cg_step (xtemp, x, d, alpha, n) ;
    cg_g (gtemp, xtemp, Com) ;
    dphi = cg_dot (gtemp, d, n) ;

    /* check if gradient is nan; if so, reduce stepsize */
    if ( dphi != dphi )
    {
        for (i = 0; i < Parm->nexpand; i++)
        {
            alpha /= rho ;
            cg_step (xtemp, x, d, alpha, n) ;
            cg_g (gtemp, xtemp, Com) ;
            dphi = cg_dot (gtemp, d, n) ;
            if ( dphi == dphi ) break ;
        }
        if ( i == Parm->nexpand )
        {
            status = -2 ;
            goto Exit ;
        }
        rho = Parm->nan_rho ;
    }
 
/*Find initial interval [a,b] such that dphia < 0, dphib >= 0,
         and phia <= phi0 + feps*fabs (phi0) */
 
    a = ZERO ;
    dphia = dphi0  ;
    ngrow = 0 ;
    nshrink = 0 ;
    while ( dphi < ZERO )
    {
        phi = cg_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) ;
            if ( phi <= fquad )
            {
                if ( Parm->PrintLevel >= 2 )
                {
                    printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n",
                            alpha, phi, fquad) ;
                }
                if ( cg_Wolfe (alpha, phi, dphi, Com) )
                {
                    status = 0 ;
                    goto Exit ;
                }
            }
        }
        if ( phi > Com->fpert )
        {
            /* 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 ;
                }
                cg_step (xtemp, x, d, alpha, n) ;
                cg_g (gtemp, xtemp, Com) ;
                dphi = cg_dot (gtemp, d, n) ;
                if ( dphi >= ZERO ) goto Secant ;
                phi = cg_f (xtemp, Com) ;
                if ( Parm->PrintLevel >= 2 )
                {
                    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 ( cg_Wolfe (alpha, phi, dphi, Com) )
                    {
                        status = 0 ;
                        goto Exit ;
                    }
                }
                if ( phi <= Com->fpert )
                {
                    a = alpha ;
                    dphia = dphi ;
                }
                else
                {
                    b = alpha ;
                }
            }
        }

/* expansion phase */

        a = alpha ;
        dphia = dphi ;
        ngrow++ ;
        if ( ngrow > Parm->nexpand )
        {
            status = 3 ;
            goto Exit ;
        }
        alpha = rho*alpha ;
        cg_step (xtemp, x, d, alpha, n) ;
        cg_g (gtemp, xtemp, Com) ;
        dphi = cg_dot (gtemp, d, n) ;
        if ( Parm->PrintLevel >= 2 )
        {
            printf ("expand,   a: %14.6e alpha: %14.6e phi: "
                     "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ;
        }
    }

Secant:
    b = alpha ;
    dphib = dphi ;
    if ( Com->QuadOK )
    {
        phi = cg_f (xtemp, Com) ;
        if ( ngrow + nshrink == 0 ) fquad = MIN (phi, Com->f0) ;
        if ( phi <= fquad )
        {
            if ( cg_Wolfe (alpha, phi, dphi, Com) )
            {
                status = 0 ;
                goto Exit ;
            }
        }
    }
    nsecant = Parm->nsecant ;
    for (iter = 1; iter <= nsecant; iter++)
    {
        if ( Parm->PrintLevel >= 2 )
        {
            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 = cg_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 >= 2 ) printf ("2nd secant\n") ;
                status = cg_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 >= 2 ) printf ("bisection\n") ;
            status = cg_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) ;
}

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

    Parm = Com->Parm ;
    if ( Parm->PrintLevel >= 1 ) printf ("Wolfe line search\n") ;
    alpha = Com->alpha ;
    phi = Com->f ;
    dphi = Com->df ;
    n = Com->n ;
    x = Com->x ;         /* current iterate */
    xtemp = Com->xtemp ; /* x + alpha*d */
    d = Com->d ;         /* current search direction */
    gtemp = Com->gtemp ; /* gradient at x + alpha*d */
    rho = Parm->rho ;
    cg_step (xtemp, x, d, alpha, n) ;
    cg_g (gtemp, xtemp, Com) ;
    dphi = cg_dot (gtemp, d, n) ;

    /* check if gradient is nan; if so, reduce stepsize */
    if ( dphi != dphi )
    {
        for (i = 0; i < Parm->nexpand; i++)
        {
            alpha /= rho ;
            cg_step (xtemp, x, d, alpha, n) ;
            cg_g (gtemp, xtemp, Com) ;
            dphi = cg_dot (gtemp, d, n) ;
            if ( dphi == dphi ) break ;
        }
        if ( i == Parm->nexpand )
        {
            status = -2 ;
            goto Exit ;
        }
        rho = Parm->nan_rho ;
    }
    dpsi = dphi - Com->wolfe_hi ;
 
/*Find initial interval [a,b] such that dphia < 0, dphib >= 0,
         and phia <= phi0 + feps*fabs (phi0) */
 
    a = ZERO ;
    dpsia = dphi0 - Com->wolfe_hi ;
    ngrow = 0 ;
    nshrink = 0 ;
    while ( dpsi < ZERO )
    {
        phi = cg_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 >= 2 )
                {
                    printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n",
                            alpha, phi, fquad) ;
                }
                if ( cg_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 ;
                }
                cg_step (xtemp, x, d, alpha, n) ;
                cg_g (gtemp, xtemp, Com) ;
                dphi = cg_dot (gtemp, d, n) ;
                dpsi = dphi - Com->wolfe_hi ;
                if ( dpsi >= ZERO ) goto Secant ;
                phi = cg_f (xtemp, Com) ;
                psi = phi - alpha*Com->wolfe_hi ;
                if ( Parm->PrintLevel >= 2 )
                {
                    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 ( cg_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 ;
        }
        alpha *= rho ;
        cg_step (xtemp, x, d, alpha, n) ;
        cg_g (gtemp, xtemp, Com) ;
        dphi = cg_dot (gtemp, d, n) ;
        dpsi = dphi - Com->wolfe_hi ;
        if ( Parm->PrintLevel >= 2 )
        {
            printf ("expand,   a: %14.6e alpha: %14.6e phi: "
                     "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ;
        }
    }

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

/* bisection iteration */

        if ( b-a >= width )
        {
            alpha = .5*(b+a) ;
            if ( Parm->PrintLevel >= 2 ) printf ("bisection\n") ;
            status = cg_updateW (&a, &dpsia, &b, &dpsib, &alpha, &phi, &dphi,
                       &dpsi, 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) ;
}

/* =========================================================================
   ==== cg_update ==========================================================
   =========================================================================
   update returns: 8 if too many iterations
                   0 if Wolfe condition is satisfied
                  -1 if interval is updated and a search is done
                  -2 if the interval updated successfully
   ========================================================================= */
int cg_update
(
    double        *a , /* left side of bracketing interval */
    double    *dphia , /* derivative at a */
    double        *b , /* right side of bracketing interval */
    double    *dphib , /* derivative at b */
    double    *alpha , /* trial step (between a and b) */
    double      *phi , /* function value at alpha (returned) */
    double     *dphi , /* function derivative at alpha (returned) */
    cg_com      *Com   /* cg com structure */
)
{
    INT n ;
    int nshrink, status ;
    double *x, *xtemp, *d, *gtemp ;
    cg_parameter *Parm ;

    Parm = Com->Parm ;
    n = Com->n ;
    x = Com->x ;         /* current iterate */
    xtemp = Com->xtemp ; /* x + alpha*d */
    d = Com->d ;         /* current search direction */
    gtemp = Com->gtemp ; /* gradient at x + alpha*d */
    cg_step (xtemp, x, d, *alpha, n) ;
    *phi = cg_fg (gtemp, xtemp, Com) ;
    *dphi = cg_dot (gtemp, d, n) ;
    if ( Parm->PrintLevel >= 2 )
    {
        printf ("update alpha: %14.6e phi: %14.6e dphi: %14.6e\n",
                 *alpha, *phi, *dphi) ;
    }
    if ( cg_Wolfe (*alpha, *phi, *dphi, Com) )
    {
        status = 0 ;
        goto Exit2 ;
    }
    status = -2 ;
    if ( *dphi >= ZERO )
    {
        *b = *alpha ;
        *dphib = *dphi ;
        goto Exit2 ;
    }
    else
    {
        if ( *phi <= Com->fpert )
        {
            *a = *alpha ;
            *dphia = *dphi ;
            goto Exit2 ;
        }
    }
    nshrink = 0 ;
    *b = *alpha ;
    while ( TRUE )
    {
        *alpha = .5*(*a + *b) ;
        nshrink++ ;
        if ( nshrink > Parm->nexpand )
        {
            status = 8 ;
            goto Exit2 ;
        }
        cg_step (xtemp, x, d, *alpha, n) ;
        *phi = cg_fg (gtemp, xtemp, Com) ;
        *dphi = cg_dot (gtemp, d, n) ;
        if ( Parm->PrintLevel >= 2 )
        {
            printf ("contract, a: %14.6e alpha: %14.6e "
                    "phi: %14.6e dphi: %14.6e\n", *a, *alpha, *phi, *dphi) ;
        }
        if ( cg_Wolfe (*alpha, *phi, *dphi, Com) )
        {
            status = 0 ;
            goto Exit2 ;
        }
        if ( *dphi >= ZERO )
        {
            *b = *alpha ;
            *dphib = *dphi ;
            goto Exit1 ;
        }
        if ( *phi <= Com->fpert )
        {
            if ( Parm->PrintLevel >= 2 )
            {
                printf ("update a: %14.6e dphia: %14.6e\n", *alpha, *dphi) ;
            }
            *a = *alpha ;
            *dphia = *dphi ;
        }
        else *b = *alpha ;
    }
Exit1:
    status = -1 ;
Exit2:
    if ( Parm->PrintLevel >= 2 )
    {
        printf ("UP a: %14.6e b: %14.6e da: %14.6e db: %14.6e status: %i\n",
                 *a, *b, *dphia, *dphib, status) ;
    }
    return (status) ;
}

/* =========================================================================
   ==== cg_updateW =========================================================
   =========================================================================
   This routine is identical to cg_update except that the function
   psi [a] = phi [a] - phi [0] - a*delta*dphi [0] is minimized instead of
   the function phi. The return int has the following meaning:
                   8 if too many iterations
                   0 if Wolfe condition is satisfied
                  -1 if interval is updated and a search is done
                  -2 if the interval updated successfully
   ========================================================================= */
int cg_updateW
(
    double        *a , /* left side of bracketing interval */
    double    *dpsia , /* derivative at a */
    double        *b , /* right side of bracketing interval */
    double    *dpsib , /* derivative at b */
    double    *alpha , /* trial step (between a and b) */
    double      *phi , /* function value at alpha (returned) */
    double     *dphi , /* derivative of phi at alpha (returned) */
    double     *dpsi , /* derivative of psi at alpha (returned) */
    cg_com      *Com   /* cg com structure */
)
{
    INT n ;
    int nshrink, status ;
    double psi, *x, *xtemp, *d, *gtemp ;
    cg_parameter *Parm ;

    Parm = Com->Parm ;
    n = Com->n ;
    x = Com->x ;         /* current iterate */
    xtemp = Com->xtemp ; /* x + alpha*d */
    d = Com->d ;         /* current search direction */
    gtemp = Com->gtemp ; /* gradient at x + alpha*d */
    cg_step (xtemp, x, d, *alpha, n) ;
    *phi = cg_fg (gtemp, xtemp, Com) ;
    psi = *phi - *alpha*Com->wolfe_hi ;
    *dphi = cg_dot (gtemp, d, n) ;
    *dpsi = *dphi - Com->wolfe_hi ;
    if ( Parm->PrintLevel >= 2 )
    {
        printf ("update alpha: %14.6e psi: %14.6e dpsi: %14.6e\n",
                 *alpha, psi, *dpsi) ;
    }
    if ( cg_Wolfe (*alpha, *phi, *dphi, Com) )
    {
        status = 0 ;
        goto Exit2 ;
    }
    status = -2 ;
    if ( *dpsi >= ZERO )
    {
        *b = *alpha ;
        *dpsib = *dpsi ;
        goto Exit2 ;
    }
    else
    {
        if ( psi <= Com->fpert )
        {
            *a = *alpha ;
            *dpsia = *dpsi ;
            goto Exit2 ;
        }
    }
    nshrink = 0 ;
    *b = *alpha ;
    while ( TRUE )
    {
        *alpha = .5*(*a + *b) ;
        nshrink++ ;
        if ( nshrink > Parm->nexpand )
        {
            status = 8 ;
            goto Exit2 ;
        }
        cg_step (xtemp, x, d, *alpha, n) ;
        *phi = cg_fg (gtemp, xtemp, Com) ;
        *dphi = cg_dot (gtemp, d, n) ;
        *dpsi = *dphi - Com->wolfe_hi ;
        psi = *phi - *alpha*Com->wolfe_hi ;
        if ( Parm->PrintLevel >= 2 )
        {
            printf ("contract, a: %14.6e alpha: %14.6e "
                    "phi: %14.6e dphi: %14.6e\n", *a, *alpha, *phi, *dphi) ;
        }
        if ( cg_Wolfe (*alpha, *phi, *dphi, Com) )
        {
            status = 0 ;
            goto Exit2 ;
        }
        if ( *dpsi >= ZERO )
        {
            *b = *alpha ;
            *dpsib = *dpsi ;
            goto Exit1 ;
        }
        if ( psi <= Com->fpert )
        {
            if ( Parm->PrintLevel >= 2 )
            {
                printf ("update a: %14.6e dpsia: %14.6e\n", *alpha, *dpsi) ;
            }
            *a = *alpha ;
            *dpsia = *dpsi ;
        }
        else *b = *alpha ;
    }
Exit1:
    status = -1 ;
Exit2:
    if ( Parm->PrintLevel >= 2 )
    {
        printf ("UP a: %14.6e b: %14.6e da: %14.6e db: %14.6e status: %i\n",
                 *a, *b, *dpsia, *dpsib, status) ;
    }
    return (status) ;
}

/* =========================================================================
   ==== cg_printParms ======================================================
   =========================================================================
   Print the contents of the cg_parameter structure
   ========================================================================= */
void cg_printParms
(
    cg_parameter  *Parm
)
{
    printf ("PARAMETERS:\n") ;
    printf ("\n") ;
    printf ("Wolfe line search parameter ..................... delta: %e\n",
             Parm->delta) ;
    printf ("Wolfe line search parameter ..................... sigma: %e\n",
             Parm->sigma) ;
    printf ("decay factor for bracketing interval ............ gamma: %e\n",
             Parm->gamma) ;
    printf ("growth factor for bracket interval ................ rho: %e\n",
             Parm->rho) ;
    printf ("growth factor for bracket interval after nan .. nan_rho: %e\n",
             Parm->nan_rho) ;
    printf ("truncation factor for cg beta ..................... eta: %e\n",
             Parm->eta) ;
    printf ("perturbation parameter for function value ......... eps: %e\n",
             Parm->eps) ;
    printf ("factor for computing average cost .............. Qdecay: %e\n",
             Parm->Qdecay) ;
    printf ("relative change in cost to stop quadstep ... QuadCufOff: %e\n",
             Parm->QuadCutOff) ;
    printf ("factor multiplying gradient in stop condition . StopFac: %e\n",
             Parm->StopFac) ;
    printf ("cost change factor, approx Wolfe transition . AWolfeFac: %e\n",
             Parm->AWolfeFac) ;
    printf ("restart cg every restart_fac*n iterations . restart_fac: %e\n",
             Parm->restart_fac) ;
    printf ("stop when cost change <= feps*|f| ................. eps: %e\n",
             Parm->eps) ;
    printf ("starting guess parameter in first iteration ...... psi0: %e\n",
             Parm->psi0) ;
    printf ("starting step in first iteration if nonzero ...... step: %e\n",
             Parm->step) ;
    printf ("factor multiply starting guess in quad step ...... psi1: %e\n",
             Parm->psi1) ;
    printf ("initial guess factor for general iteration ....... psi2: %e\n",
             Parm->psi2) ;
    printf ("max iterations is n*maxit_fac ............... maxit_fac: %e\n",
             Parm->maxit_fac) ;
    printf ("max expansions in line search ................. nexpand: %i\n",
             Parm->nexpand) ;
    printf ("max secant iterations in line search .......... nsecant: %i\n",
             Parm->nsecant) ;
    printf ("print level (0 = none, 2 = maximum) ........ PrintLevel: %i\n",
             Parm->PrintLevel) ;
    printf ("Logical parameters:\n") ;
    if ( Parm->PertRule )
        printf ("    Error estimate for function value is eps\n") ;
    else
        printf ("    Error estimate for function value is eps*Ck\n") ;
    if ( Parm->QuadStep )
        printf ("    Use quadratic interpolation step\n") ;
    else
        printf ("    No quadratic interpolation step\n") ;
    if ( Parm->PrintFinal )
        printf ("    Print final cost and statistics\n") ;
    else
        printf ("    Do not print final cost and statistics\n") ;
    if ( Parm->PrintParms )
        printf ("    Print the parameter structure\n") ;
    else
        printf ("    Do not print parameter structure\n") ;
    if ( Parm->AWolfe)
        printf ("    Approximate Wolfe line search\n") ;
    else
        printf ("    Wolfe line search") ;
        if ( Parm->AWolfeFac > 0. )
            printf (" ... switching to approximate Wolfe\n") ;
        else
            printf ("\n") ;
    if ( Parm->StopRule )
        printf ("    Stopping condition uses initial grad tolerance\n") ;
    else
        printf ("    Stopping condition weighted by absolute cost\n") ;
    if ( Parm->debug)
        printf ("    Check for decay of cost, debugger is on\n") ;
    else
        printf ("    Do not check for decay of cost, debugger is off\n") ;
}

/*
Version 1.2 Change:
    The variable dpsi needs to be included in the argument list for
    subroutine cg_updateW (update of a Wolfe line search)

Version 2.0 Changes:
    The user interface was redesigned. The parameters no longer need to
    be read from a file. For compatibility with earlier versions of the
    code, we include the routine cg_readParms to read parameters.
    In the simplest case, the user can use NULL for the
    parameter argument of cg_descent, and the code sets the default
    parameter values. If the user wishes to modify the parameters, call
    cg_default in the main program to initialize a cg_parameter
    structure. Individual elements of the structure could be modified.
    The header file cg_user.h contains the structures and prototypes
    that the user may need to reference or modify, while cg_descent.h
    contains header elements that only cg_descent will access.  Note
    that the arguments of cg_descent have changed.

Version 3.0 Changes:
    Major overhaul
*/