========================================================================= */int asa_descent /*  return:                      -5 (ginorm < tau2*pgnorm without hitting boundary)                      -4 (ginorm >=tau2*pgnorm, many x components hit boundary)                      -3 (ginorm >=tau2*pgnorm, one x component hits boundary)                      -2 (function value became nan)                      -1 (starting function value is nan)                       0 (convergence tolerance satisfied)                       1 (change in func <= feps*|f|)                       2 (total iterations exceeded maxit)                       3 (slope always negative in line search)                       4 (number secant iterations exceed nsecant)                       5 (search direction not a descent direction)                       6 (line search fails in initial interval)                       7 (line search fails during bisection)                       8 (line search fails during interval update)                       9 (debugger is on and the function value increases)*/(    asa_com *Com){    int     i, iter, j, maxit, n, n5, nfree, nf, ng, nrestart,            status ;    double  delta2, eta_sq, Qk, Ck, pgnorm, ginorm,            f, ftemp, gnorm, xnorm, gnorm2, dnorm2, denom,            t, t1, t2, t3, t4, t5, dphi, dphi0, alpha, talpha,            xj, gj, xg, xp, sts, sty, sk,            yk, ykyk, ykgk, dkyk, yk1, yk2, yk3, yk4, yk5, beta,            *x, *d, *g, *xtemp, *gtemp, *lo, *hi, *pg ;    asacg_parm *Parm ;    asa_parm *asaParm ;/* initialization */    x = Com->x ;    lo = Com->lo ;    hi = Com->hi ;    n = Com->n ;    d = Com->d ;    g = Com->g ;    xtemp = Com->xtemp ;    gtemp = Com->gtemp ;    pg = Com->pg ;    nfree = Com->nfree ;    nf = Com->nf ;    ng = Com->ng ;    pgnorm = Com->pgnorm ;    ginorm = Com->ginorm ;    Parm = Com->cgParm ;    asaParm = Com->asaParm ;    if ( Parm->PrintLevel >= 1 )    {        printf ("Dimension in CG, nfree = %i\n", nfree) ;    }    /* the conjugate gradient algorithm is restarted every nrestart iteration */    nrestart = (INT32) (((double) nfree)*Parm->restart_fac) ;    /* abort when number of iterations reaches maxit in one pass through cg */    if ( Parm->maxit_fac == INF ) Com->cgmaxit = maxit = INT_INF ;    else Com->cgmaxit = maxit = (INT32) (((double) n)*Parm->maxit_fac) ;    n5 = nfree % 5 ;    f = Com->f ;    Ck = ZERO ;    Qk = ZERO ;/* initial function and gradient evaluations, initial direction */    Com->f0 = f + f ;    xnorm = asa_max (x, nfree) ;    gnorm = ZERO ;    gnorm2 = ZERO ;    for (i = 0; i < n5; i++)    {        t = g [i] ;        d [i] = -t ;        gnorm2 += t*t ;        if ( gnorm < fabs (t) ) gnorm = fabs (t) ;    }    for (; i < nfree;)    {        t1 = g [i] ;        d [i] = -t1 ;        if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ;        i++ ;        t2 = g [i] ;        d [i] = -t2 ;        if ( gnorm < fabs (t2) ) gnorm = fabs (t2) ;        i++ ;        t3 = g [i] ;        d [i] = -t3 ;        if ( gnorm < fabs (t3) ) gnorm = fabs (t3) ;        i++ ;        t4 = g [i] ;        d [i] = -t4 ;        if ( gnorm < fabs (t4) ) gnorm = fabs (t4) ;        i++ ;        t5 = g [i] ;        d [i] = -t5 ;        if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ;        i++ ;        gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ;    }    /* check that starting function value is nan */    if ( f != f )    {        status = -1 ;        goto Exit ;    }    if ( Parm->PrintLevel >= 2 )    {        printf ("iter: %5i f = %14.6e pgnorm = %14.6e ginorm = %14.6e\n\n",          (int) 0, f, pgnorm, ginorm) ;    }    dphi0 = -gnorm2 ;    delta2 = 2*Parm->delta - ONE ;    eta_sq = Parm->eta*Parm->eta ;    alpha = Parm->step ;    if ( alpha == ZERO )    {        alpha = Parm->psi0*xnorm/gnorm ;        if ( xnorm == ZERO )        {            if ( f != ZERO ) alpha = Parm->psi0*fabs (f)/gnorm2 ;            else             alpha = ONE ;        }    }/*  start the conjugate gradient iteration    alpha starts as old step, ends as final step for current iteration    f is function value for alpha = 0    Com->QuadOK = TRUE means that a quadratic step was taken */    for (iter = 1; iter <= maxit; iter++)    {        Com->QuadOK = FALSE ;        alpha = Parm->psi2*alpha ;        asa_maxstep (x, d, Com) ;        if ( Parm->QuadStep )        {            if ( f != ZERO ) t = fabs ((f-Com->f0)/f) ;            else             t = ONE ;            if ( t > Parm->QuadCutOff )       /* take provisional step talpha */            {                talpha = Parm->psi1*alpha ;                if ( Com->minstep >= asaParm->parm1*talpha )                {                    talpha = MIN (talpha, asaParm->parm2*Com->minstep) ;                    asa_step (xtemp, x, d, talpha, nfree) ;                    /*provisional function value*/                    ftemp = asa_f (xtemp, Com) ;                    /* check if function value is nan */                    if ( ftemp != ftemp ) /* reduce stepsize */                    {                        for (i = 0; i < Parm->nexpand; i++)                        {                            talpha /= Parm->rho ;                            asa_step (xtemp, x, d, talpha, nfree) ;                            ftemp = asa_f (xtemp, Com) ;                            if ( ftemp == ftemp ) break ;                        }                        if ( i == Parm->nexpand )                        {                            status = -2 ;                            goto Exit ;                        }                    }                    if ( ftemp < f )              /* check if quadstep > 0 */                    {                       denom = TWO*(((ftemp-f)/talpha)-dphi0) ;                       if ( denom > ZERO )    /* try a quadratic fit step */                       {                           Com->QuadOK = TRUE ;                           alpha = -dphi0*talpha/denom ;                       }                    }                }            }        }        Com->f0 = f ;                          /* f0 saved as prior value */        if ( Parm->PrintLevel >= 3 )        {            printf ("minstep =%14.6e, maxstep =%14.6e \n",                     Com->minstep, Com->maxstep) ;            printf ("QuadOK: %2i initial a: %14.6e f0: %14.6e dphi0: %14.6e\n",                    Com->QuadOK, alpha, Com->f0, dphi0) ;            if ( (alpha > Com->minstep) && Com->QuadOK )            {                printf("Quadratic step > minstep to boundary\n") ;            }        }        /* parameters in Wolfe, approximate Wolfe conditions, and in update */        Qk = Parm->Qdecay*Qk + ONE ;        Ck = Ck + (fabs (f) - Ck)/Qk ;        /* average cost magnitude */        if ( Parm->PertRule ) Com->fpert = f + Parm->eps*Ck ;        else                  Com->fpert = f + Parm->eps ;        Com->wolfe_hi = Parm->delta*dphi0 ;        Com->wolfe_lo = Parm->sigma*dphi0 ;        Com->awolfe_hi = delta2*dphi0 ;        Com->alpha = alpha ;/* either double prior step or quadratic fit step */        Com->f = f ;        if ( Com->AWolfe )                  /* approximate Wolfe line search*/        {            if ( Parm->PrintLevel >= 3 )            {                printf ("Perform approximate Wolfe line search\n") ;            }            status = asa_line (dphi0, Com) ;        }        else                                  /* ordinary Wolfe line search */        {            if ( Parm->PrintLevel >= 3 )            {                 printf ("Perform ordinary Wolfe line search\n") ;            }            status = asa_lineW (dphi0, Com) ;        }        /* if ordinary Wolfe line search fails, possibly try approximate           Wolfe line search*/        if ( (status > 0) && !Com->AWolfe && (Parm->AWolfeFac > ZERO) )        {            Com->AWolfe = TRUE ;            if ( Parm->PrintLevel >= 3 )            {                printf ("Ordinary Wolfe line search fails, "                        "try approximate Wolfe line search\n") ;            }            status = asa_line (dphi0, Com) ;        }        alpha = Com->alpha ;        f = Com->f ;        dphi = Com->df ;        if ( (status > 0) || (status == -1) || (status == -2) ) goto Exit ;        /*Test for convergence to within machine epsilon          [set feps to zero to remove this test] */        if ( (-alpha*dphi0 <= Parm->feps*fabs (f)) && (status == 0) )        {            status = 1 ;            goto Exit ;        }        /* compute beta, yk2, gnorm, gnorm2, dnorm2, update x and g */        if ( iter % nrestart != 0 )        {            ginorm = ZERO ;            pgnorm = ZERO ;            for (j = 0; j < nfree; j++)            {                xj = xtemp [j] ;                gj = gtemp [j] ;                xg = xj - gj ;                if      ( xg > hi [j] ) xp = hi [j] - xj ;                else if ( xg < lo [j] ) xp = xj - lo [j] ;                else                    xp = fabs (gj) ;                pgnorm = MAX (pgnorm, xp) ;                ginorm = MAX (ginorm, fabs (gj)) ;            }            for (; j < n; j++)            {                xj = x [j] ;                gj = gtemp [j] ;                xg = xj - gj ;                if      ( xg > hi [j] ) xp = hi [j] - xj ;                else if ( xg < lo [j] ) xp = xj - lo [j] ;                else                    xp = fabs (gj) ;                pgnorm = MAX (pgnorm, xp) ;            }            if ( asa_tol (pgnorm, Com) )            {                status = 0 ;                for (j = 0; j < nfree; j++)                {                    xj = xtemp [j] ;                    x [j] = xj ;                    gj = gtemp [j] ;                    xg = xj - gj ;                    g [j] = gj ;                    if      ( xg > hi [j] ) pg [j] = hi [j] - xj ;                    else if ( xg < lo [j] ) pg [j] = lo [j] - xj ;                    else                    pg [j] = -gj ;                }                for (; j < n; j++)                {                    xj = x [j] ;                    gj = gtemp [j] ;                    xg = xj - gj ;                    g [j] = gj ;                    if      ( xg > hi [j] ) pg [j] = hi [j] - xj ;                    else if ( xg < lo [j] ) pg [j] = lo [j] - xj ;                    else                    pg [j] = -gj ;                }                goto Exit1 ;            }            if ( ginorm < pgnorm*Com->tau2 ) status = -5 ;            if ( status < -2 )            {                sts = ZERO ;                sty = ZERO ;                for (j = 0; j < nfree; j++)                {                    t = xtemp[j] ;                    sk = t - x [j] ;                    x [j] = t ;                    sts += sk*sk ;                    t = gtemp [j] ;                    sty += sk*(t-g [j]) ;                    g [j] = t ;                }                Com->sts = sts ;                Com->sty = sty ;                goto Exit ;            }            asa_copy (x, xtemp, nfree) ;            dnorm2 = ZERO ;            for (j = 0; j < n5; j++) dnorm2 = dnorm2 + d [j]*d [j] ;            for (; j < nfree; j += 5)            {                dnorm2 = dnorm2 + d [j]*d [j] + d [j+1]*d [j+1]                                              + d [j+2]*d [j+2]                                              + d [j+3]*d [j+3]                                              + d [j+4]*d [j+4] ;            }            ykyk = ZERO ;            ykgk = ZERO ;            for (j = 0; j < n5; j++)            {                t = gtemp [j] ;                yk = t - g [j] ;                ykyk += yk*yk ;                ykgk += yk*t ;                g [j] = t ;            }            for (j = n5; j < nfree; )            {                t1 = gtemp [j] ;                yk1 = t1 - g [j] ;                g [j] = t1 ;                j++ ;                t2 = gtemp [j] ;                yk2 = t2 - g [j] ;                g [j] = t2 ;                j++ ;                t3 = gtemp [j] ;                yk3 = t3 - g [j] ;                g [j] = t3 ;                j++ ;                t4 = gtemp [j] ;                yk4 = t4 - g [j] ;                g [j] = t4 ;                j++ ;                t5 = gtemp [j] ;                yk5 = t5 - g [j] ;                g [j] = t5 ;                j++ ;                ykyk += yk1*yk1 + yk2*yk2 + yk3*yk3 + yk4*yk4 + yk5*yk5 ;                ykgk += yk1*t1 + yk2*t2 + yk3*t3 + yk4*t4 + yk5*t5 ;            }            dkyk = dphi - dphi0 ;