/* =========================================================================   ============================== ASA_CG ===================================   =========================================================================       ________________________________________________________________      | A Conjugate Gradient (cg_descent) based Active Set Algorithm   |      |                                                                |      |                  Version 1.0  (May 18, 2008)                   |      |                                                                |      |        William W. Hager    and      Hongchao Zhang             |      |        hager@math.ufl.edu         hzhang@math.ufl.edu          |      |  Department of Mathematics      Department of Mathematics      |      |     University of Florida       Louisiana State University     |      |  Gainesville, Florida 32611      Baton Rouge, Louisiana        |      |     352-392-0281 x 244                                         |      |                                                                |      |      Copyright by William W. Hager and Hongchao Zhang          |      |                                                                |      |          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 "asa_user.h"#include <iostream>#include "asa_cg.h"using namespace std;static int iii = 1;int asa_cg /*  return:                      -2 (function value became nan in cg)                      -1 (starting function value is nan in cg)                       0 (convergence tolerance satisfied)                       1 (change in func <= feps*|f| in cg)                       2 (cg iterations in all passes or                          in one pass exceeded their limit)                       3 (slope always negative in line search in cg)                       4 (number secant iterations exceed nsecant in cg)                       5 (search direction not a descent direction in cg)                       6 (line search fails in initial interval in cg)                       7 (line search fails during bisection in cg)                       8 (line search fails during interval update in cg)                       9 (debugger is on and the function value increases in cg)                      10 (out of memory)                      11 (cbb iterations in all passes or                          in one pass exceeded their limit)                      12 (line search failed in cbb iteration)                      13 (search direction in cbb is not descent direction)                      14 (function value became nan in cbb) */(    double            *x, /* input: starting guess, output: the solution */    double           *lo, /* lower bounds */    double           *hi, /* upper bounds */    INT32                n, /* problem dimension */    asa_stat       *Stat, /* structure with statistics (can be NULL) */    asacg_parm    *CParm, /* user parameters, NULL = use default parameters */    asa_parm      *AParm, /* user parameters, NULL = use default parameters */    double      grad_tol, /* |Proj (x_k - g_k) - x_k|_inf <= grad_tol */    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),                          NULL = compute value & gradient using value & grad */    double        *Work  /* either work array of size 7n + memory (m) or NULL */){    int gp, ident, j, nfree, status, *ifree ;    INT32 cbb_totit, cg_totit ;    double alpha, gj, pert_lo, pert_hi, t, tl, th, gnorm, ginorm, pgnorm, xnorm,           xj, xg, xp, *work, *d, *g, *xtemp, *gtemp, *pg ;    asacg_parm *cgParm, cgParmStruc ;    asa_parm *asaParm, asaParmStruc ;    asa_com Com ;/* initialize the parameters */    if ( CParm == NULL )    {        cgParm = &cgParmStruc ;        asa_cg_default (cgParm) ;    }    else cgParm = CParm ;    if ( cgParm->PrintParms ) asa_printcgParms (cgParm) ;    if ( AParm == NULL )    {        asaParm = &asaParmStruc ;        asa_default (asaParm) ;    }    else asaParm = AParm ;    if ( asaParm->PrintParms ) asa_printParms (asaParm) ;    /* abort after maxit iterations of cbb in one pass */    if ( asaParm->maxit_fac == INF ) Com.pgmaxit = INT_INF ;    else Com.pgmaxit = (INT32) (((double) n)*asaParm->maxit_fac) ;    /* abort after totit iterations of cbb in all passes */    if ( asaParm->totit_fac == INF ) cbb_totit = INT_INF ;    else cbb_totit = (INT32) (((double) n)*asaParm->totit_fac) ;    /* abort after maxfunc function evaluation in one pass of cbb */    if ( asaParm->maxfunc_fac == INF ) Com.pgmaxfunc = INT_INF ;    else Com.pgmaxfunc = (INT32) (((double) n)*asaParm->maxfunc_fac) ;    /* abort after totit iterations of cg in all passes */    if ( cgParm->totit_fac == INF ) cg_totit = INT_INF ;    else cg_totit = (INT32) (((double) n)*cgParm->totit_fac) ;    pert_lo = asaParm->pert_lo ;    pert_hi = asaParm->pert_hi ;    Com.tau1 = asaParm->tau1 ;    Com.tau2 = asaParm->tau2 ;    Com.cgParm = cgParm ;    Com.asaParm = asaParm ;    Com.x = x ;    Com.n = n ;             /* problem dimension */    Com.n5 = n % 5 ;    Com.nf = (INT32) 0 ;      /* number of function evaluations */    Com.ng = (INT32) 0 ;      /* number of gradient evaluations */    Com.cbbiter = (INT32) 0 ; /* number of cbb iterations evaluations */    Com.cgiter = (INT32) 0 ;  /* number of cg iterations */    Com.AWolfe = cgParm->AWolfe ; /* do not touch user's AWolfe */    Com.AArmijo = asaParm->AArmijo ; /* do not touch user's AArmijo */    Com.value = value ;    Com.grad = grad ;    Com.valgrad = valgrad ;    Com.DimReduce = FALSE ;    ifree = Com.ifree = (int*)malloc (n*sizeof (int)) ;    if ( Work == NULL ) work = (double*)malloc ((5*n+asaParm->m)*sizeof (double)) ;    else                work = Work ;    if ( work == NULL )    {        printf ("Insufficient memory for specified problem dimension %e\n",                 (double) n) ;        status = 10 ;        return (status) ;    }    d = Com.d = work ;    g = Com.g = d+n ;    xtemp = Com.xtemp = g+n ;    gtemp = Com.gtemp = xtemp+n ;    pg = Com.pg = gtemp+n ;    Com.lastfvalues = pg+n ; /* size asaParm->m */    Com.lo = lo ;    Com.hi = hi ;    Com.cbbiter = 0 ;    Com.cbbfunc = 0 ;    Com.cbbgrad = 0 ;    Com.cgiter = 0 ;    Com.cgfunc = 0 ;    Com.cggrad = 0 ;    ident = FALSE ;    xnorm = ZERO ;    for (j = 0; j < n; j++)    {        t = x [j] ;        if      ( t > hi [j] ) t = hi [j] ;        else if ( t < lo [j] ) t = lo [j] ;        x [j] = t ;        if ( xnorm < fabs (t) ) xnorm = fabs (t) ;    }    Com.f = asa_fg (g, x, &Com) ;    pgnorm = ZERO ;    gnorm = ZERO ;    for (j = 0; j < n; j++)    {        xj = x [j] ;        gj = g [j] ;        xg = xj - gj ;        if      ( xg > hi [j] ) xp = hi [j] - xj ;        else if ( xg < lo [j] ) xp = lo [j] - xj ;        else                    xp = -gj ;        pg [j] = xp ;        pgnorm = MAX (pgnorm, fabs (xp)) ;        gnorm = MAX (gnorm, fabs (gj)) ;    }    if ( asaParm->StopRule ) Com.tol = MAX (pgnorm*asaParm->StopFac, grad_tol) ;    else                     Com.tol = grad_tol ;    Com.pgnorm = Com.pgnorm_start = pgnorm ;    if ( asa_tol (pgnorm, &Com) )    {        status = 0 ;        goto Exit ;    }    if ( xnorm != ZERO ) Com.alpha = alpha = xnorm/gnorm ;    else                 Com.alpha = alpha = ONE/gnorm ;    /* compute gradient norm for inactive variables */    ginorm = ZERO ;    nfree = 0 ;    gp = FALSE ;    for (j = 0; j < n; j++)    {        xj = x [j] ;        tl = lo [j] ;        th = hi [j] ;        gj = g [j] ;        xg = xj - alpha*gj ;        if      ( (xg >= th) && (th-xj > pert_hi) ) gp = TRUE ;        else if ( (xg <= tl) && (xj-tl > pert_lo) ) gp = TRUE ;        if ( (xj-tl > pert_lo) && (th - xj > pert_hi) )        {            ginorm = MAX (ginorm, fabs (gj)) ;            ifree [nfree] = j ;            nfree++ ;        }    }    Com.ginorm = ginorm ;    Com.nfree = nfree ;    if ( asaParm->PrintLevel >= 1 )    {        printf ("\ninitial f = %14.6e pgnorm = %14.6e ginorm = %14.6e\n",                 Com.f, pgnorm, ginorm) ;        printf ("            nfree = %i xnorm = %14.6e gp = %i\n",                 nfree, xnorm, gp) ;    }    if ( (ginorm < Com.tau1*pgnorm) || gp || asaParm->GradProjOnly )    {        Com.cbbfunc = 1 ;        Com.cbbgrad = 1 ;        goto Grad_proj ;    }    else    {        Com.cgfunc = 1 ;        Com.cggrad = 1 ;        goto CG_descent ;    }    Grad_proj:    if ( asaParm->PrintLevel >= 1 ) printf ("\nGradProj:\n") ;    Com.DimReduce = FALSE ;    status = asa_grad_proj(&Com) ;    if ( asaParm->PrintLevel >= 1 )    {        printf ("exit Grad_proj\n") ;    }    if ( Com.cbbiter >= cbb_totit ) status = 11 ;    if ( status >= 0 ) goto Exit ;    /* extract free variable */    nfree = 0 ;    for (j = 0; j < n; j++)    {        xj = x [j] ;        if ( (xj-lo [j] > pert_lo) && (hi [j] - xj > pert_hi) )        {            ifree [nfree] = j ;            nfree++ ;        }    }    Com.nfree = nfree ;    CG_descent:    if ( nfree != n )    {       asa_shrink_all (&Com) ;       asa_copy (xtemp+nfree, x+nfree, n-nfree) ;       Com.DimReduce = TRUE ;    }    else Com.DimReduce = FALSE ;    if ( asaParm->PrintLevel >= 1 ) printf ("\nCG:\n") ;    status = asa_descent (&Com) ;    if ( asaParm->PrintLevel >= 1 )    {        printf ("exit the CG subroutine\n") ;    }    if ( Com.DimReduce ) asa_expand_all (&Com) ;    if ( Com.cgiter >= cg_totit ) status = 2 ;    if ( status >= -2 ) goto Exit ;    /* ginorm < tau2* pgnorm without hitting boundary */    if ( status == -5 )    {        Com.alpha = asa_init_bbstep (&Com) ;        goto Grad_proj ;    }    /* ginorm >= tau2* pgnorm and many components of x hit boundary  */    else if ( status == -4 )    {        ginorm = ZERO ;        nfree = 0 ;        for (j = 0 ; j < n; j++)        {            xj = x [j] ;            if ( (xj-lo [j] > pert_lo) && (hi [j] - xj > pert_hi) )            {                t = fabs (g [j]) ;                ginorm = MAX (ginorm, t) ;                ifree [nfree] = j ;                nfree++ ;            }        }        Com.nfree = nfree ;        Com.ginorm = ginorm ;        if ( ginorm >= Com.tau1*Com.pgnorm ) goto CG_descent ;        else        {           if ( asaParm->PrintLevel >= 1 ) printf ("ginorm < tau1* pgnorm\n") ;           Com.alpha = asa_init_bbstep (&Com) ;           goto Grad_proj ;        }    }    /* ginorm >= tau2* pgnorm and only one component of x hits boundary */    else if ( status == -3 )    {        if ( pgnorm < asaParm->pgdecay*MAX (Com.pgnorm_start, ONE) )        {            ident = asa_identify (x, g, Com.pgnorm, &Com) ;        }        if ( ident )        {            ident = FALSE ;            ginorm = ZERO ;            nfree = 0 ;            for (j = 0 ; j < n; j++)            {                xj = x [j] ;                if ( (xj-lo [j] > pert_lo) && (hi [j] - xj > pert_hi) )                {                    t = fabs (g [j]) ;                    ginorm = MAX (ginorm, t) ;                    ifree [nfree] = j ;                    nfree++ ;                }            }            Com.nfree = nfree ;            Com.ginorm = ginorm ;            if ( ginorm >= Com.tau1*Com.pgnorm ) goto CG_descent ;            else            {               if ( asaParm->PrintLevel >= 1 )                   printf ("ginorm < tau1* pgnorm\n" ) ;               Com.alpha = asa_init_bbstep (&Com) ;               goto Grad_proj ;            }        }        else        {            Com.alpha = asa_init_bbstep (&Com) ;            goto Grad_proj ;        }    }    Exit:    if ( (asaParm->PrintFinal) || (asaParm->PrintLevel >= 1) )    {        const char mess1 [] = "Possible causes of this error message:" ;        const char mess2 [] = "   - your tolerance may be too strict: "                              "grad_tol = " ;        const char mess4 [] = "   - your gradient routine has an error" ;        const char mess5 [] = "   - the parameter epsilon in "                              "asa_descent_c.parm is too small" ;        printf ("\nFinal convergence status = %d\n", status);        if ( status == -2 )        {            printf ("Function value became nan at cg iteration %10.0e\n",                     (double) Com.cgiter) ;        }        else if ( status == -1 )        {            printf ("Function value of starting point is nan at "                     "cg iteration %10.0f\n", (double) Com.cgiter) ;        }        else if ( status == 0 )        {