C     =================================================================
C     Module: Augmented Lagrangian subroutines
C     =================================================================

C     Last update of any of the component of this module: 
C
C     September 13, 2006.

C     ******************************************************************
C     ******************************************************************

      subroutine easyalabb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml,
     +meq,alin,acol,aval,annz,b,epsgap,epsfeas,epsopt,maxoutit,iprint,
     +ncomp,macheps,bignum,f,snorm,nalpsupn,ldh,H,outiter,totnnodes,
     +inform,wi1,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,wd11,wd12)

      implicit none

C     This subroutine aims to simplify the usage of algbetra.

C     SCALAR ARGUMENTS
      integer annz,inform,iprint,ldh,m,maxoutit,meq,ml,n,ncomp,outiter,
     +        totnnodes
      double precision bignum,epsgap,epsfeas,epsopt,f,macheps,nalpsupn,
     +        snorm

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz),wi1(ml)
      logical equatn(m),linear(m),linvar(n)
      double precision aval(annz),b(ml),l(n),lambda(m),H(ldh,n),rho(m),
     +        u(n),wd1(m),wd2(m),wd3(n),wd4(m),wd5(m),wd6(n),wd7(n),
     +        wd8(n),wd9(n),wd10(n),wd11(n),wd12(n),x(n)

C     LOCAL SCALARS
      logical rhoauto
      character pptype * 3
      integer maxitncp
      double precision epsgaini,epsgafin,epsgacon,rhomult,rhofrac,
     +        lammin,lammax

C     SET AUGMENTED-LAGRANGIAN PARAMETERS

      lammin    = - 1.0d+20         
      lammax    =   1.0d+20         

      rhoauto   =    .true. 
      rhomult   =    10.0d0
      rhofrac   =     0.5d0
      pptype    =     'ONE'

      epsgaini  =   1.0d-01
      epsgacon  =   1.0d-01
      epsgafin  =    epsgap

      maxitncp  =         9

C     CALL ALABB

      call alabb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,alin,
     +acol,aval,annz,b,lammin,lammax,rhoauto,rhomult,rhofrac,pptype,
     +epsgaini,epsgafin,epsgacon,epsfeas,epsopt,maxitncp,maxoutit,
     +iprint,ncomp,macheps,bignum,f,wd1,wd2,snorm,wd3,nalpsupn,ldh,H,
     +outiter,totnnodes,inform,wd4,wd5,wi1,wd6,wd7,wd8,wd9,wd10,wd11,
     +wd12)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine alabb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,
     +alin,acol,aval,annz,b,lammin,lammax,rhoauto,rhomult,rhofrac,
     +pptype,epsgaini,epsgafin,epsgacon,epsfeas,epsopt,maxitncp,
     +maxoutit,iprint,ncomp,macheps,bignum,f,c,sigma,snorm,nal,nalpsupn,
     +ldh,H,outiter,totnnodes,inform,lambar,sigpre,wi1,wd1,wd2,wd3,wd4,
     +wd5,wd6,wd7)

      implicit none

C     Solves the general-constrained minimization problem
C
C     min f(x)
C
C     subject to
C
C             c_j(x)  = 0, j in E,
C             c_j(x) <= 0, j in I,
C             l <= x <= u,
C
C     where E is the set of indices of the equality constraints, I is
C     the set of indices of the inequality constraints, and there are
C     n variables and m constraints, using the method of multipliers 
C     described in:
C 
C     R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, 
C     "On Augmented Lagrangian methods with general lower-level 
C     constraints", Technical Report MCDO-050304, Department of Applied 
C     Mathematics, UNICAMP, Brazil, 2005.

C     SCALAR ARGUMENTS
      logical rhoauto
      character * 3 pptype
      integer annz,inform,iprint,ldh,m,maxitncp,maxoutit,meq,ml,n,ncomp,
     +        outiter,totnnodes
      double precision bignum,epsfeas,epsgaini,epsgafin,epsgacon,epsopt,
     +        f,lammax,lammin,macheps,nalpsupn,rhofrac,rhomult,snorm

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz),wi1(ml)
      logical equatn(m),linear(m),linvar(n)
      double precision aval(annz),b(ml),c(m),H(ldh,n),l(n),lambar(m),
     +        lambda(m),nal(n),rho(m),sigma(m),sigpre(m),u(n),wd1(n),
     +        wd2(n),wd3(n),wd4(n),wd5(n),wd6(n),wd7(n),x(n)

C     Parameters of the subroutine:
C     =============================
C
C     On Entry:
C     =========
C
C     N integer: Number of variables.
C     ---------
C
C     X double precision x(n): Initial estimation of the solution.
C     -----------------------
C
C     L double precision l(n): Lower bounds for x.
C     -----------------------
C
C     U double precision u(n): Upper bounds for x.
C     -----------------------
C
C     M integer
C     ---------
C
C     Number of constraints (equalities plus inequalities) without 
C     considering the bound constraints.
C
C     LAMBDA double precision lambda(m)
C     ---------------------------------
C
C     Initial estimation of the Lagrange multipliers.
C
C     RHO double precision rho(m): Initial penalty parameters.
C     ---------------------------
C
C     EQUATN logical equatn(m)
C     ------------------------
C
C     Logical array that, for each constraint, indicates whether the 
C     constraint is an equality constraint (TRUE) or an inequality 
C     constraint (FALSE).
C
C     LINEAR logical linear(m)
C     ------------------------
C
C     Logical array that, for each constraint, indicates whether the 
C     constraint is a linear constraint (TRUE) or a nonlinear constraint 
C     (FALSE).
C
C     LAMMIN double precision: See below.
C     -----------------------
C
C     LAMMAX double precision
C     -----------------------
C
C     The Lagrange multipliers approximations used for the definition of 
C     the subproblems must be inside the interval [lammin,lammax].
C
C     RHOMULT double precision
C     ------------------------
C
C     When increased, the penalty parameters will be multiplied by 
C     rhomult.
C
C     RHOFRAC double precision
C     ------------------------
C
C     The penalty parameters will be, basically, increased when the 
C     infeasibility is larger than a fraction of the previous 
C     infeasibility. This fraction is rhofrac.
C
C     PPTYPE character * 3
C     --------------------
C
C     Indicates the type of penalty parameters that will be used. 
C     pptype = 'ONE' means just a unique penalty parameter for all the 
C     constraints while pptype = 'DIF' means a penalty parameter per 
C     constraint.
C
C     EPSFEAS double precision 
C     ------------------------
C
C     Required precision for feasibility and complementarity.
C
C     EPSOPT double precision
C     ----------------------- 
C
C     Required precision for optimality.
C
C     MAXITNCP integer 
C     ----------------
C
C     Maximum number of outer iterations withoutprogress in feasibility.
C
C     MAXOUTIT integer
C     ----------------
C
C     Maximum number of Augmented Lagrangian (outer) iterations.
C     (Ignored in the unconstrained and bound-constrained cases.)
C
C     IPRINT integer
C     --------------
C                
C     Controls the ammount of information of the output according to the 
C     following convention:
C
C     0 means no output at all.
C
C     1 means information at each outer iteration but without any 
C       information of how the subproblems are being solved.
C
C     2 means the same as 1 plus information of each inner iteration.
C
C     3 means the same as 2 plus information of the line searches and 
C       the calculation of the truncated Newton direction (using CG) of 
C       each inner iteration.
C
C     In all cases, an output file named solution.txt with the final 
C     point, Lagrange mutipliers and penalty parameters will be 
C     generated. Moreover, the same output of the screen will be saved 
C     in a file named alabb.out.
C
C     NCOMP integer
C     -------------
C
C     Every time a vector is printed, just its first ncomp component 
C     will be displayed.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     BIGNUM double precision: A big number like 1.0d+99.
C     -----------------------
C
C     LAMBAR and SIGPRE double precision lambar(m),sigpre(m)
C     ------------------------------------------------------
C
C     M-dimensional double precision working vectors.
C
C     WI1 integer wi1(n): N-dimensional integer working vector.
C     ------------------
C
C     WD1 ... WD14 double precision wd1(n) ... wd14(n)
C     ------------------------------------------------
C
C     N-dimensional double precision working vectors.
C              
C     On Return:
C     ==========
C
C     X double precision x(n): Final estimation of the solution.
C     -----------------------
C
C     LAMBDA double precision lambda(m)
C     ---------------------------------
C
C     Final estimation of the Lagrange multipliers.
C
C     RHO double precision rho(m): Final penalty parameters.
C     ---------------------------
C
C     F double precision: Objective functional value at the solution.
C     ------------------
C
C     C double precision c(m): Constraints at the solution.
C     -----------------------
C
C     SIGMA double precision sigma(m)
C     -------------------------------
C
C     Complementarity measurement at the solution: c(i) for equality 
C     constraints and max(c(i), - lambar(i) / rho(i)) for inequalities.
C
C     SNORM double precision: Sup-norm of the constraints.
C     ----------------------
C
C     NAL double precision nal(n)
C     ---------------------------
C
C     Gradient of the Augmented Lagrangian at the solution.
C
C     NALPSUPN double precision
C     -------------------------
C 
C     Sup-norm of the continuous projected gradient of the Augmented 
C     Lagrangian at the solution
C
C     OUTITER integer: Number of outer iterations.
C     ---------------
C 
C     INFORM integer
C     --------------
C
C     This output parameter tells what happened in this subroutine, 
C     according to the following conventions:
C 
C       0 means convergence with feasibility, optimality and 
C         complementarity.
C
C       1 means that it was achieved the maximum number of outer 
C         iterations (maxoutit).
C
C       4 means that the algorithm stopped by ``lack of feasibility 
C         progress'', i.e., the current point is infeasible (the 
C         constraints violation is larger than the tolerance epsfeas) 
C         and the constraint violation has not even simple decrease 
C         during maxitncp consecutive iterations. In this case, the 
C         problem may be infeasible.
C
C     -90 means that subroutine evalf   retuned an error flag.
C
C     -91 means that subroutine evalc   retuned an error flag.
C
C     -92 means that subroutine evalg   retuned an error flag.
C
C     -93 means that subroutine evaljac retuned an error flag.
C
C     -94 means that subroutine evalh   retuned an error flag.
C 
C     -95 means that subroutine evalhc  retuned an error flag.

C     LOCAL SCALARS
      integer i,icrit,mprint,nprint,nnodes,rhoincr
      double precision al,snormprev,epsgak,gap,nalpi,rhoini,sumc

C     SET UP SOME CONSTANTS

C     just for printing
      nprint   = min0(n,ncomp)
      mprint   = min0(m,ncomp)

C     INTIALIZATION

      icrit     = 0
      outiter   = 0
      totnnodes = 0

      gap = bignum

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      call evalobjc(n,x,f,m,c,inform)
      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

C     SET PENALTY PARAMETERS:

C     AUGMENTED LAGRANGIAN FUNCTION IS OF THE FORM
C
C         AL = F(X) + \SUM_i P(LAMBDA_i,C_i(X),RHO),
C
C     WHERE
C
C         P(LAMBDA,Y,RHO) = Y ( LAMBDA + 0.5 RHO Y ), 
C
C     IF C_i(X) IS AN EQUALITY CONSTRAINT OR LAMBDA + RHO Y > 0, AND
C
C         P(LAMBDA,Y,RHO) = - 0.5 LAMBDA^2 / RHO,
C
C     OTHERWISE.
C
C     ASSUMING THAT LAMBDA_i = 0 FOR ALL i, IT IS CLEAR THAT 
C
C         P(LAMBDA_i,C_i(X),RHO) = 0.5 RHO C_i(X)^2
C
C     AND THAT THE VALUE OF RHO THAT BALANCES F(X) AND 
C     \SUM_i P(LAMBDA_i,C_i(X),RHO) IS GIVEN BY
C
C     RHO = F(X) / ( 0.5 \SUM_i C_i(X)^2 ),
C
C     WHERE THE SUM IS MADE OVER THE VIOLATED CONSTRAINTS.

      if ( rhoauto ) then

          sumc = 0.0d0
          do i = 1,m
              if ( equatn(i) .or. c(i) .gt. 0.0d0 ) then
                  sumc = sumc + 0.5d0 * c(i) ** 2
              end if
          end do

          if ( sumc .eq. 0.0d0 ) then
              rhoini = 10.0d0
          else
              rhoini = 10.0d0 * max( 1.0d0, abs( f ) ) / sumc
              rhoini = max( 1.0d-06, min( rhoini, 10.0d0 ) )
          end if 

          do i = 1,m
              rho(i) = rhoini
          end do

      end if

C     COMPUTE LAGRANGE MULTIPLIERS

      do i = 1,m
          lambar(i) = max( lammin, min( lammax, lambda(i) ) )
      end do

C     COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS
 
      snormprev = bignum

      snorm = 0.0d0
      do i = 1,m
          if ( equatn(i) ) then
              sigma(i) = c(i)
          else
              sigma(i) = max( c(i), - lambar(i) / rho(i) )
          end if
          snorm = max( snorm, abs( sigma(i) ) )
      end do

C     COMPUTE CONTINUOUS PROJECTED GRADIENT NORM

      call evalnal(n,x,m,lambar,rho,equatn,linear,nal,inform)
      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

      nalpsupn = 0.0d0
      do i = 1,n
          nalpi    = max( l(i), min( x(i) - nal(i), u(i) ) ) - x(i)
          nalpsupn = max( nalpsupn, abs( nalpi ) )
      end do

C     PRINT INITIAL INFORMATION

      if  ( iprint .ge. 1 ) then
          write(*,  900) n,m
          write(*,  910) nprint,(l(i),i=1,nprint)
          write(*,  920) nprint,(u(i),i=1,nprint)

          write(10, 900) n,m
          write(10, 910) nprint,(l(i),i=1,nprint)
          write(10, 920) nprint,(u(i),i=1,nprint)
      end if

C     MAIN LOOP

 100  continue

C     PRINT ITERATION INFORMATION

      if ( iprint .ge. 1 ) then
          write(*,  970) outiter
          write(*,  980) nprint,(x(i),i=1,nprint)
          write(*,  990) mprint,(lambda(i),i=1,mprint)
          write(*, 1000) mprint,(rho(i),i=1,mprint)
          write(*, 1020) f,snorm,nalpsupn

          write(10, 970) outiter
          write(10, 980) nprint,(x(i),i=1,nprint)
          write(10, 990) mprint,(lambda(i),i=1,mprint)
          write(10,1000) mprint,(rho(i),i=1,mprint)
          write(10,1020) f,snorm,nalpsupn
      end if

C     SAVING INTERMEDIATE DATA FOR CRASH REPORT

      open(20,file='alabb-tabline.out')
      write(20,3000) 0.0d0,-1,n,m,ml,outiter,totnnodes,f,snorm,nalpsupn
      close(20)

C     TEST STOPPING CRITERIA

      if ( gap .le. epsgafin .and. 
     +     snorm .le. epsfeas .and. nalpsupn .le. epsopt ) then

        ! THE POINT IS FEASIBLE AND OPTIMAL

          inform = 0

          if ( iprint .ge. 1 ) then
              write(*, 1060) inform,epsfeas,epsfeas,epsopt
              write(10,1060) inform,epsfeas,epsfeas,epsopt
          end if

          go to 500
      end if

      if ( outiter .ge. maxoutit ) then
        ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP
          inform = 1

          if ( iprint .ge. 1 ) then
              write(*, 1070) inform,maxoutit
              write(10,1070) inform,maxoutit
          end if

          go to 500
      end if
 
C     TEST IF THE INFEASIBILITY DECREASED SUFFICIENTLY

      if ( snorm .gt. epsfeas .and. snorm .ge. snormprev ) then

          icrit = icrit + 1

          if ( icrit .ge. maxitncp ) then
            ! LACK OF FEASIBILITY PROGRESS, STOP
              inform = 4

              if ( iprint .ge. 1 ) then
                  write(*, 1100) inform,maxitncp
                  write(10,1100) inform,maxitncp
              end if

              go to 500
          end if

      else
          icrit = 0
      end if

C     DO AN OUTER ITERATION

      outiter = outiter + 1      

C     OPTIMALITY REQUERIMENT FOR THE SUBPROBLEM

      if ( outiter .eq. 1 ) then
          epsgak = epsgaini
      else
          epsgak = max( epsgafin, epsgak * epsgacon )
          if ( epsgak .le. 2.0d0 * epsgafin ) epsgak = epsgafin
      end if

C     CALL THE INNER-SOLVER

      call abb(n,x,l,u,m,lambar,rho,equatn,linear,linvar,ml,meq,alin,
     +acol,aval,annz,b,epsgak,epsopt,iprint,ncomp,macheps,bignum,al,nal,
     +nalpsupn,gap,ldh,H,nnodes,inform,wi1,wd1,wd2,wd3,wd4,wd5,wd6,wd7)

      totnnodes = totnnodes + nnodes

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

C     COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS

      call evalobjc(n,x,f,m,c,inform)
      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

C     UPDATE LAGRANGE MULTIPLIERS APPROXIMATION 

      do i = 1,m
          call evaldpdy(c(i),rho(i),lambar(i),equatn(i),lambda(i))
          lambar(i) = max( lammin, min( lammax, lambda(i) ) )
      end do

C     COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS
 
      snormprev = snorm

      snorm = 0.0d0
      do i = 1,m
          sigpre(i) = sigma(i)

          if ( equatn(i) ) then
              sigma(i) = c(i)
          else
              sigma(i) = max( c(i), - lambar(i) / rho(i) )
          end if

          snorm = max( snorm, abs( sigma(i) ) )
      end do

C     UPDATE PENALTY PARAMETERS

      if ( pptype .eq. 'ONE' ) then

          if ( outiter .gt. 1 .and. 
     +         snorm .gt. max( epsfeas, rhofrac * snormprev ) ) then

              do i = 1,m
                  rho(i) = rhomult * rho(i)
              end do

              if ( iprint .ge. 1 ) then
                  write(*, 1030) rhofrac * snormprev,rhomult
                  write(10,1030) rhofrac * snormprev,rhomult
              end if

          else

              if ( iprint .ge. 1 ) then
                  write(*, 1050)
                  write(10,1050)
              end if

          end if

      else ! if ( pptype .eq. 'DIF' ) then

          rhoincr = 0
          do i = 1,m
              if ( outiter .gt. 1 .and. abs(sigma(i)) .gt. 
     +        max ( epsfeas, rhofrac * abs(sigpre(i)) ) ) then
                  rho(i)  = rhomult * rho(i)
                  rhoincr = rhoincr + 1
              end if
          end do

          if ( iprint .ge. 1 ) then

              if ( rhoincr .ne. 0 ) then
                  write(*, 1040) rhoincr,m,rhomult
                  write(10,1040) rhoincr,m,rhomult
              else
                  write(*, 1050)
                  write(10,1050)
              end if

          end if

      end if

C     ITERATE

      go to 100

C     STOP

 500  continue

C     NON-EXECUTABLE STATEMENTS 

  900 format( /,' Entry to ALABB.',/,
     +        /,' Number of variables: ',I7,
     +          ' Number of constraints: ',I7)
  910 format( /,' Lower bounds (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  920 format( /,' Upper bounds (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  970 format( /,' ALABB outer iteration: ',I7)
  980 format( /,' Current point (first ',I7, ' components): ',
     +        /,6(1X,1PD11.4))
  990 format( /,' Updated Lagrange multipliers (first ',I7,
     +          ' components): ',
     +        /,6(1X,1PD11.4))
 1000 format( /,' Updated penalty parameters (first ',I7,
     +          ' components): ',
     +        /,6(1X,1PD11.4))
 1020 format( /,' Objective function value',38X,' = ',1PD11.4,
     +        /,' Sup-norm of the constraints',35X,' = ',1PD11.4,
     +        /,' Sup-norm of the projected gradient of the',
     +          ' Lagrangian',10X,' = ',1PD11.4)
 1030 format( /,' The desired infeasibility was not achieved (',
     +            1PD11.4,'). The penalty',
     +        /,' parameter will be increased multiplying by rhofrac',
     +          ' = ',1PD11.4,'.')
 1040 format( /,' The desired infeasibility was not achieved in ',I7,
     +          ' over ',I7,
     +        /,' constraints.',
     +        /,' Penalty parameters will be increased multiplying',
     +          ' by rhofrac = ',1PD11.4,'.')
 1050 format( /,' Desired feasibility improvement was achieved,',
     +        /,' so, penalty parameter(s) will not be modified.')
 1060 format( /,' Flag of ALABB = ',I2,
     +          ' (Convergence with Sup-norm of the Lagrangian',
     +          ' projected',/,' gradient smaller than ',1PD11.4,',',
     +          ' constraints Sup-norm smaller than',
     +        /,' ',1PD11.4,', and largest complementarity violation',
     +          ' smaller than ',1PD11.4,')')
 1070 format( /,' Flag of ALABB = ',I2,
     +          ' (It was exceeded the maximum allowed number of outer',
     +        /,' iterations (maxoutit =',I7,').)')
 1100 format( /,' Flag of ALABB = ',I2,
     +          ' (Constraint violations have not decreased',
     +          ' substantially',/,' over',I2,' outer iterations.',
     +          ' Problem possibly infeasible.)')
 2000 format( /,' Flag of ALABB = ',I3,' Fatal Error',/,
     +        /,' The following codes means: ',/,
     +        /,' -90 : error in evalf   subroutine', 
     +        /,' -91 : error in evalc   subroutine', 
     +        /,' -92 : error in evalg   subroutine', 
     +        /,' -93 : error in evaljac subroutine', 
     +        /,' -94 : error in evalh   subroutine', 
     +        /,' -95 : error in evalhc  subroutine',/) 
 3000 format(F8.2,1X,I3,1X,I6,1X,I6,1X,I6,1X,I3,1X,I7,1X,1P,D24.16,1X,
     +       1P,D7.1,1X,1P,D7.1)

      end

