printf ("Convergence tolerance for gradient satisfied\n") ;        }        else if ( status == 1 )        {            printf ("Terminating in cg since change in function value "                    "<= feps*|f|\n") ;        }        else if ( status == 2 )        {            printf ("Number of iterations exceed specified limits "                    "for cg routine\n") ;            printf ("Iterations: %10.0f maxit: %10.0f totit: %10.0f\n",                    (double) Com.cgiter, (double) Com.cgmaxit,                    (double) cg_totit) ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;        }        else if ( status == 3 )        {            printf ("Slope always negative in cg 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 in cg, too many secant steps\n") ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;        }        else if ( status == 5 )        {            printf ("Search direction not a descent direction in cg\n") ;        }        else if ( status == 6 ) /* line search fails */        {            printf ("Line search fails in cg iteration\n") ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;            printf ("%s\n", mess4) ;            printf ("%s\n", mess5) ;        }        else if ( status == 7 ) /* line search fails */        {            printf ("Line search fails in cg iteration\n") ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;        }        else if ( status == 8 ) /* line search fails */        {            printf ("Line search fails in cg iteration\n") ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;            printf ("%s\n", mess4) ;            printf ("%s\n", mess5) ;        }        else if ( status == 9 )        {            printf ("Debugger is on, function value does not improve in cg\n") ;            printf ("new value: %25.16e old value: %25.16e\n",                Com.f_debug, Com.f0) ;        }        else if ( status == 10 )        {            printf ("Insufficient memory\n") ;        }        else if ( status == 11 )        {            printf ("Number of iterations or function evaluation exceed "                          "specified limits for cbb routine\n") ;            printf ("Iterations: %i maxit: %i totit: %i\n",                     Com.cbbiter, Com.pgmaxit, cbb_totit) ;            printf ("Total function evaluations: %i maxfunc: %i\n",                     Com.nf, Com.pgmaxfunc);        }        if ( status == 12 ) /* line search fails in cbb iteration */        {            printf ("Line search fails in cbb iteration\n") ;            printf ("%s\n", mess1) ;            printf ("%s %e\n", mess2, Com.tol) ;            printf ("%s\n", mess4) ;        }        if ( status == 13 )        {            printf ("Search direction not descent direction in "                    "asa_grad_proj\n") ;            printf ("directional derivative: %e\n", Com.gtd) ;        }        if ( status == 14 )        {             printf ("At cbb iteration %i function value became nan\n",                      Com.cbbiter) ;        }        printf ("projected gradient max norm: %13.6e\n", Com.pgnorm) ;        printf ("function value:              %13.6e\n", Com.f) ;        printf ("\nTotal cg  iterations:           %10.0f\n",                (double) Com.cgiter) ;        printf ("Total cg  function evaluations: %10.0f\n",                (double) Com.cgfunc) ;        printf ("Total cg  gradient evaluations: %10.0f\n",                (double) Com.cggrad) ;        printf ("Total cbb iterations:           %10.0f\n",                (double) Com.cbbiter) ;        printf ("Total cbb function evaluations: %10.0f\n",                (double) Com.cbbfunc) ;        printf ("Total cbb gradient evaluations: %10.0f\n",                    (double) Com.cbbgrad) ;        printf ("------------------------------------------\n") ;        printf ("Total function evaluations:     %10.0f\n",                (double) Com.nf) ;        printf ("Total gradient evaluations:     %10.0f\n",                (double) Com.ng) ;        printf ("==========================================\n\n") ;    }    free (ifree) ;    if ( Work == NULL ) free (work) ;    if ( Stat != NULL )    {        Stat->f = Com.f ;        Stat->pgnorm = Com.pgnorm ;        Stat->cgiter = Com.cgiter ;        Stat->cgfunc = Com.cgfunc ;        Stat->cggrad = Com.cggrad ;        Stat->cbbiter = Com.cbbiter ;        Stat->cbbfunc = Com.cbbfunc ;        Stat->cbbgrad = Com.cbbgrad ;    }    return (status) ;}/* =========================================================================   === asa_default ======================================================   =========================================================================   Set default parameter values for the ASA routine. The CG default   parameter values are set by asa_cg_default.  If the parameter argument of   asa_descent is NULL, this routine is called by asa_cg automatically.   If the user wishes to set parameter values, then the asa_parameter structure   should be allocated in the main program. The user could call asa_default   to initialize the structure, and then individual elements in the structure   could be changed, before passing the structure to asa_cg.   =========================================================================*/void asa_default(    asa_parm *Parm){    double eps, t ;    /* T => print final statistics       F => no printout of statistics */    Parm->PrintFinal = TRUE ;    /* Level 0  = no printing), ... , Level 4 = maximum printing */    Parm->PrintLevel = 0 ;    /* T => print parameters values       F => do not display parmeter values */    Parm->PrintParms = FALSE ;    /* T => use approximate nonmonotone Armijo line search       F => use ordinary nonmonotone Armijo line search, switch to            approximate Armijo when |f_r-f| < AArmijoFac*|min (f_r, f_{max})| */    Parm->AArmijo = FALSE ;    Parm->AArmijoFac = 1.e-8 ;    /* Stop Rules (these override the corresponding cg parameters):       T => ||proj_grad||_infty <= max(grad_tol,initial ||grad||_infty*StopFact)       F => ||proj_grad||_infty <= grad_tol*(1 + |f_k|) */    Parm->StopRule = TRUE ;    Parm->StopFac = 0.e-12 ;    /* T => estimated error in function value = eps*|min (f_r, f_{max}) |       F => estimated error in function value = eps */    Parm->PertRule = TRUE ;    Parm->eps = 1.e-6 ;    /* T => only use gradient projection algorithm       F => let algorithm decide between grad_proj and cg_descent */    Parm->GradProjOnly = FALSE ;    /* abort cbb after maxit_fac*n iterations in one pass through cbb */    Parm->maxit_fac = INF ;    /* abort cbb after totit_fac*n iterations in all passes through cbb */    Parm->totit_fac = INF ;    /* abort cbb iteration after maxfunc_fac*n function evaluations */    Parm->maxfunc_fac = INF ;    /* perturbation in bounds based on machine epsilon, which we now compute */    eps = ONE ;    t = ONE ;    while ( t > 0 )    {        eps /= TWO ;        t = ONE + eps ;        t -= ONE ;    }    eps *= 2 ;                   /* machine epsilon */    Parm->pert_lo = 1.e3*eps ;   /* perturbation of lower bounds */    Parm->pert_hi = 1.e3*eps ;   /* perturbation of upper bounds */    /* search for non nan function value by shrinking search interval       at most nshrink times */    Parm->nshrink = (int) 50 ;    /* factor by which interval shrinks when searching for non nan value */    Parm->nan_fac = 2.e-1 ;    /* update fr if fmin was not improved after L iterations */    Parm->L = 3 ;    /* fmax = max (f_{k-i}, i = 0, 1, ..., min (k, m-1) ) */    Parm->m = 8 ;    /* update fr if initial stepsize was accepted in previous P iterations */    Parm->P = 40 ;    /* CBB cycle length */    Parm->nm = 4 ;    /* Reinitialize BB stepsize, if (s^t y)/(||s|| ||y||) >= gamma       and ||s|| <= min (parm3*|f_k+1|/||g_k+1||_infty, 1) */    Parm->gamma = 0.975e0 ;    /* update reference value fr if (fr-fmin)/(fc-fmin) > gamma1 */    Parm->gamma1 = (double) Parm->m / (double) Parm->L ;    /* update fr if (fr-f)/(fmax-f) > gamma2, np > P, and fmax > f */    Parm->gamma2 = (double) Parm->P / (double) Parm->m ;    /* terminate Armijo line search when       phi(alpha) <= phi_r + alpha * delta * phi'(0) where phi_r = fr or fcomp*/    Parm->delta = 1.0e-4 ;   /* Armijo line search parameter */    /* stepsize s in the line search must satisfy lmin <= s <= lmax */    Parm->lmin = 1.0e-20 ;    Parm->lmax = 1.0e+20 ;    /* attempt a quadratic interpolation step in cg_descent if the       provisional stepsize times parm1 <= stepsize to boundary */    Parm->parm1 = 1.e-1 ;    /* if quadratic interpolation step is attempted, the provisional step       is at most parm2*stepsize to boundary */    Parm->parm2 = 9.e-1 ;    /* used in the the criterion of reinitializing the BB stepsize */    Parm->parm3 = 1.e-1 ;    /* maximum number of previous BB steps used when s^t y <= ZERO */    Parm->parm4 = 6 ;    /* if ginorm < tau1*pgnorm, continue gradient projection steps  */    Parm->tau1 = 1.e-1 ;    /* decay factor for tau1 */    Parm->tau1_decay = 5.e-1 ;    /* ginorm < tau2*pgnorm implies subproblem solved in cgdescent */    Parm->tau2 = 1.e-1 ;    /* decay factor for tau2 */    Parm->tau2_decay = 5.e-1 ;    /* if pgnorm < pgdecay*MAX (pgnorm0, ONE), check the undecided index set                                pgnorm0 = pgnorm at starting point */    Parm->pgdecay = 1.e-4 ;    /* backtracking decay factor in the Armijo line search */    Parm->armijo_decay = 5.e-1 ;    /* use quadratic interpolation to compute Armijo step if it       lies in the interval [.1 alpha, .9 alpha] */    Parm->armijo0 = 1.e-1 ;    Parm->armijo1 = 9.e-1 ;}/* =========================================================================   === asa_cg_default ======================================================   =========================================================================   Set default conjugate gradient parameter values. If the parameter argument   of asa_cg is NULL, this routine is called by asa_cg automatically.   If the user wishes to set parameter values, then the asa_parameter structure   should be allocated in the main program. The user could call asa_cg_default   to initialize the structure, and then individual elements in the structure   could be changed, before passing the structure to asa_cg.   =========================================================================*/void asa_cg_default(    asacg_parm   *Parm){    /* Level 0 = no printing, ... , Level 4 = 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 ;    /* T => estimated error in function value is eps*Ck,       F => estimated error in function value is eps */    Parm->PertRule = TRUE ;    Parm->eps = 1.e-4 ;    /* 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 ;    /* 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 ;    /* if step is nonzero, it is the initial step of the initial line search */    Parm->step = ZERO ;    /* abort cg after maxit_fac*n iterations in one pass */    Parm->maxit_fac = INF ;    /* abort cg after totit_fac*n iterations in all passes */    Parm->totit_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. ;}/* =========================================================================   === asa_descent =========================================================   =========================================================================   cg_descent conjugate gradient algorithm with modifications to handle the   bound constraints.