#include <math.h>#include <limits.h>#include <float.h>#include <string.h>#include <ctype.h>#include <stdlib.h>#include <stdio.h>#define ZERO ((double) 0)#define ONE  ((double) 1)#define TWO  ((double) 2)#define MAX(a,b) (((a) > (b)) ? (a) : (b))#define MIN(a,b) (((a) < (b)) ? (a) : (b))typedef struct asa_com_struct /* common variables */{    /* parameters computed by the code */    double          *lo ; /* lower bounds */    double          *hi ; /* lower bounds */    double           *x ; /* current estimate for solution */    double           *d ; /* current search direction */    double           *g ; /* current gradient */    double       *xtemp ; /* x + alpha*d */    double       *gtemp ; /* gradient at x + alpha*d */    double          *pg ; /* projected gradient */    double *lastfvalues ; /* previous function values */    double      minstep ; /* smallest step to reach a bound */    double      maxstep ; /* step which makes all bounds active */    int           nfree ; /* number of free variables */    int          *ifree ; /* indices of free variables */    int               n ; /* problem dimension, saved for reference */    int              n5 ; /* n % 5 */    INT32              nf ; /* total number of function evaluations */    INT32              ng ; /* total number of gradient evaluations */    INT32         cbbiter ; /* total number of cbb iterations */    INT32          cgiter ; /* total number of cg iterations */    INT32         cbbfunc ; /* total cbb function evaluations */    INT32         cbbgrad ; /* total cbb gradient evaluations */    INT32          cgfunc ; /* total cg function evaluations */    INT32          cggrad ; /* total cg gradient evaluations */    INT32         cgmaxit ; /* maximum number of iterations in one pass of cg */    INT32         pgmaxit ; /* max iterations in one pass of gradient projection*/    INT32       pgmaxfunc ; /* max function evals/pass of gradient projection */    double        alpha ; /* stepsize along search direction */    double            f ; /* function value for step alpha */    double      f_debug ; /* function value at time of debug failure */    double           df ; /* function derivative for step alpha */    double        fpert ; /* perturbation is eps*Ck if PertRule is T */    double           f0 ; /* old function value */    double           Ck ; /* average cost as given by the rule:                             Qk = Qdecay*Qk + 1, Ck += (fabs (f) - Ck)/Qk */    double     wolfe_hi ; /* upper bound for slope in Wolfe test */    double     wolfe_lo ; /* lower bound for slope in Wolfe test */    double    awolfe_hi ; /* upper bound for slope, approximate Wolfe test */    int          QuadOK ; /* T (quadratic step successful) */    int          AWolfe ; /* F (use Wolfe line search)                             T (use approximate Wolfe line search)                             do not change user's AWolfe, this value can be                             changed based on AWolfeFac */    int       DimReduce ;  /*T (compressed problem, nfree < n)                             F (work in full space, nfree = n)*/    int         AArmijo ; /* F (use nonmonotone Armijo line search)                             T (use approximate nonmonotome Armijo line search)                             do not change user's AArmijo, this value can be                             changed based on AArmijoFac */    double          sts ; /* ||s||^2 */    double          gtd ; /* g'd */    double          sty ; /* s'y */    double      pert_lo ; /* perturbation of lower bounds */    double      pert_hi ; /* perturbation of upper bounds */    double         tau1 ; /* if ginorm < tau1*pgnorm, continue gp steps  */    double         tau2 ; /* ginorm < tau2*pgnorm => subproblem solved in cgdescent */    double pgnorm_start ; /* ||Proj (x_0 - g_0) - x_0||_infty */    double          tol ; /* convergence tolerance */    double       pgnorm ; /* project gradient norm */    double       ginorm ; /* norm of inactive gradient components */    asacg_parm  *cgParm ; /* cg user parameters */    asa_parm   *asaParm ; /* asa user parameters */    double        (*value) (double *, INT32) ; /* f = value (x, n) */    void           (*grad) (double *, double *, INT32) ; /* grad (g, x, n) */    double      (*valgrad) (double *, double *, INT32) ; /* f = valgrad (g,x,n)*/} asa_com ;/* prototypes */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 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 */) ;int asa_tol(    double      pgnorm, /* projected gradient sup-norm */    asa_com       *Com) ;double asa_dot /* dot product between x and y */(    double *x , /* first vector */    double *y , /* second vector */    int     n   /* length of vectors */) ;void asa_copy  /* Copy vector x into vector y */(    double *y,    double *x,    int     n) ;void asa_saxpy  /* Compute z = y + ax */(    double *z,    double *y,    double *x,    double  a,    int     n) ;double asa_max /* Return max {fabs (x [j]) : 1 <= j < n}*/(    double *x,    int     n) ;void asa_step /* Compute xtemp = x + alpha d */(    double *xtemp , /*output vector */    double     *x , /* initial vector */    double     *d , /* search direction */    double  alpha , /* stepsize */    INT32         n   /* length of the vectors */) ;int asa_line(    double       dphi0, /* function derivative at starting point (alpha = 0) */    asa_com       *Com  /* cg com structure */) ;int asa_lineW(    double       dphi0 , /* function derivative at starting point (alpha = 0) */    asa_com       *Com   /* cg com structure */) ;int asa_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) */    asa_com       *Com   /* cg com structure */) ;int asa_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) */    asa_com       *Com   /* cg com structure */) ;void asa_project(    double  *xnew,    double     *x,    double     *d,    double  alpha,    asa_com  *Com   /* cg com structure */) ;void asa_maxstep(    double       *x, /* current iterate */    double       *d, /* direction */    asa_com    *Com) ;int asa_grad_proj /*return:                      -1 (proceed to cg routine)                       0 (convergence tolerance satisfied)                      10 (number of iterations exceeds limit)                      11 (line search fails) */(    asa_com *Com) ;double asa_f(    double    *x,    asa_com *Com) ;void asa_g(    double    *g,    double    *x,    asa_com *Com) ;double asa_fg(    double    *g,    double    *x,    asa_com *Com) ;int asa_identify(   double     *x,   double     *g,   double pgnorm,   asa_com  *Com) ;double asa_init_bbstep(    asa_com *Com) ;void asa_expandx(    double    *x,    asa_com *Com) ;void asa_shrinkx(    double    *x,    asa_com *Com) ;void asa_shrinkxg(    double    *x,    double    *g,    asa_com *Com) ;void asa_expand_all(    asa_com *Com) ;void asa_shrink_all(    asa_com *Com) ;void asa_printcgParms(    asacg_parm  *Parm) ;void asa_printParms(    asa_parm  *Parm) ;