C     =================================================================
C     File: alabb.f
C     =================================================================

C     =================================================================
C     Module: Auxiliary subroutines
C     =================================================================

C     Last update of any of the component of this module: 
C
C     March 27, 2006.

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

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

      implicit none

C     This subroutine: 
C
C     (1) Computes some machine-dependent constants. 
C
C     (2) Process the optionally input specification file alabb.dat. 
C
C     (3) Check analytic derivatives if desired by the user.
C
C     (4) Open and close the output file alabb.out. 
C
C     (5) Computes the CPU elapsed time used by the solver.
C
C     (6) Calls the solver. 
C
C     (7) Write an output file called solution.txt with the solution.

C     SCALAR ARGUMENTS
      logical checkder
      integer annz,inform,iprint,ldh,m,maxoutit,meq,ml,n,ncomp,outiter,
     +        totnnodes
      double precision epsgap,epsfeas,epsopt,f,gap,nalpsupn,snorm
      real time

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz),wi1(n),wi2(ml)
      logical equatn(m),linear(m),linvar(n)
      double precision aval(annz),b(ml),H(ldh,n),l(n),lambda(m),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     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      integer i
      double precision bignum,macheps

C     DECLARATIONS RELATED TO THE TIME MEASUREMENT USING DTIME
      real dtime
      external dtime

      real dum(2)
      data dum/0.0,0.0/

C     EXTERNAL FUNCTIONS
      double precision d1mach

C     COMPUTE MACHINE-DEPENDENT CONSTANTS

      bignum = 1.0d+99

      macheps = d1mach(4)

C     OPEN THE OUTPUT FILE

      open(10,file='alabb.out')

C     SET SOME SOLVER ARGUMENTS USING THE SPECIFICATION FILE

      call fparam(checkder,epsgap,epsfeas,epsopt,maxoutit,iprint,ncomp)

C     TEST DERIVATIVES

      if ( checkder ) then
          call checkd(m,n,l,u,wi1,wd3,wd6,wd7,wd8,wd9,macheps,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
      end if

C     CALL ALABB IF THE PROBLEM HAS CONSTRAINTS. OTHERWISE, CALL ABB

      time = dtime(dum)

      if ( m .gt. 0 ) then

C         CALL THE AUGMENTED LAGRANGIAN SOLVER ALABB

          call 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,wi2,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,
     +    wd11,wd12)

      else

C         CALL THE BOX-CONSTRAINTS GLOBAL SOLVER ABB

          call abb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,
     +    alin,acol,aval,annz,b,epsgap,epsopt,iprint,ncomp,macheps,
     +    bignum,f,wd3,nalpsupn,gap,ldh,H,totnnodes,inform,wi2,wd4,wd5,
     +    wd6,wd7,wd8,wd9,wd10,wd11,wd12)

          outiter =     0
          snorm   = 0.0d0

      end if

      time = dtime(dum)
      time = dum(1)

C     CLOSE THE OUTPUT FILE

      close(10)

C     SAVE THE SOLUTION

      open(10,file='solution.txt')

C     Solution point
      write(10,9000)
      do i = 1,n
          write(10,9010) i,x(i)
      end do

C     Lagrange multipliers and penalty parameters
      if ( m .gt. 0 ) then
          write(10,9020)
          do i = 1,m
              write(10,9030) i,lambda(i),rho(i)
          end do
      end if

      close(10)

C     WRITE STATISTICS

      open(10,file='alabb-tabline.out')

      write(10,9050) time,inform,n,m,ml,outiter,totnnodes,f,snorm,
     +nalpsupn

      close(10)

C     RETURN

 500  continue

C     NON-EXECUTABLE STATEMENTS

 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',/) 
 9000 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)')
 9010 format(  I7,1P,D24.16)
 9020 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS AND ',
     +         'PENALTY PARAMETERS: ',//,2X,'INDEX',11X,
     +         'LAMBDA(INDEX)',15X,'RHO(INDEX)')
 9030 format(I7,1P,D24.16,1X,1P,D24.16)
 9050 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

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

      subroutine fparam(checkder,epsgap,epsfeas,epsopt,maxoutit,iprint,
     +ncomp)

      implicit none

C     This subroutine set some alabb arguments related to stopping 
C     criteria and output. The setting values are taken fron file 
C     alabb.dat. Nothing is done if the file does not exist. File 
C     alabb.dat is an alternative to run ALABB several times 
C     with different arguments without having to compile it again.

C     SCALAR ARGUMENTS
      logical checkder
      integer iprint,maxoutit,ncomp
      double precision epsfeas,epsgap,epsopt

C     PARAMETERS
      integer nwords
      parameter ( nwords = 7 )

C     DATA BLOCKS
      character * 1 lower(26),upper(26)
      data lower /'a','b','c','d','e','f','g','h','i','j','k','l','m',
     +            'n','o','p','q','r','s','t','u','v','w','x','y','z'/
      data upper /'A','B','C','D','E','F','G','H','I','J','K','L','M',
     +            'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/

      character * 10 dictionary(nwords)
      data dictionary(1) /'CHECK-DERI'/
      data dictionary(2) /'GAP-TOLERA'/
      data dictionary(3) /'FEASIBILIT'/
      data dictionary(4) /'OPTIMALITY'/
      data dictionary(5) /'MAX-OUTER-'/
      data dictionary(6) /'OUTPUT-DET'/
      data dictionary(7) /'NCOMP-ARRA'/

      character * 32 description(nwords)
      data description(1) /'CHECK-DERIVATIVES               '/
      data description(2) /'GAP-TOLERANCE                   '/
      data description(3) /'FEASIBILITY-TOLERANCE           '/
      data description(4) /'OPTIMALITY-TOLERANCE            '/
      data description(5) /'MAX-OUTER-ITERATIONS            '/
      data description(6) /'OUTPUT-DETAIL                   '/
      data description(7) /'NCOMP-ARRAY                     '/

      character * 1 addinfo(nwords)
      data addinfo(1) /' '/
      data addinfo(2) /'D'/
      data addinfo(3) /'D'/
      data addinfo(4) /'D'/
      data addinfo(5) /'I'/
      data addinfo(6) /'I'/
      data addinfo(7) /'I'/

C     LOCAL SCALARS
      integer i,ifirst,ikey,ilast,inum,j
      double precision dnum

C     LOCAL ARRAYS
      character * 80 line
      character * 10 keyword

C     OPENING THE SPECIFICATION FILE
      open(20,err=300,file='alabb.dat',status='old')

C     MAIN LOOP
      write(*, 9005)
      write(10,9005)

 100  continue

C     READING LINES
      read(20,fmt=1000,err=400,end=200) line

C     PROCESS LINES

C     Find first character
      i = 1
 110  if ( i .le. 80 .and. line(i:i) .eq. ' ' ) then
          i = i + 1
          go to 110
      end if

C     Skip blank lines
      if ( i .gt. 80 ) then
          go to 100
      end if

      ifirst = i

C     Skip comments
      if ( line(ifirst:ifirst) .eq. '*' .or. 
     +     line(ifirst:ifirst) .eq. '#' ) then
          go to 100
      end if

C     Find the end of the keyword
      i = ifirst + 1
 120  if ( i .le. 80 .and. line(i:i) .ne. ' ' ) then
          i = i + 1
          go to 120
      end if

      ilast = i - 1

C     Obtain the first 10 characters and convert to upper-case
      keyword = '          '
      do i = 1,min( 10, ilast - ifirst + 1 )
          keyword(i:i) = line(ifirst+i-1:ifirst+i-1)
          do j = 1,26
              if ( keyword(i:i) .eq. lower(j) ) then
                  keyword(i:i) = upper(j)
              end if
          end do
      end do

C     Look up the keyword in the dictionary
      i = 1
 130  if ( i .le. nwords .and. keyword .ne. dictionary(i) ) then
          i = i + 1
          go to 130
      end if

C     Ignore unknown keywords
      if ( i .gt. nwords ) then
          write(*, 9020) line(ifirst:ilast)
          write(10,9020) line(ifirst:ilast)
          go to 100
      end if

      ikey = i

C     Read additional information if needed
      if ( addinfo(ikey) .ne. ' ' ) then

C         Skip blanks
          i = ilast + 1
 140      if ( i .le. 80 .and. line(i:i) .eq. ' ' ) then
              i = i + 1
              go to 140
          end if

C         Ignore keywords without the required information
          if ( i .gt. 80 ) then
              write(*, 9030) description(ikey)
              write(10,9030) description(ikey)
              go to 100
          end if 

C         Read additional information
          if ( addinfo(ikey) .eq. 'I' ) then
              read(unit=line(i:80),fmt=2000) inum

          else if ( addinfo(ikey) .eq. 'D' ) then
              read(unit=line(i:80),fmt=3000) dnum
          end if

      end if

C     Process keyword
      if ( addinfo(ikey) .eq. ' ' ) then
          write(*, 9040) description(ikey)
          write(10,9040) description(ikey)
      else if ( addinfo(ikey) .eq. 'I' ) then
          write(*, 9041) description(ikey),inum
          write(10,9041) description(ikey),inum
      else if ( addinfo(ikey) .eq. 'D' ) then
          write(*, 9042) description(ikey),dnum
          write(10,9042) description(ikey),dnum
      end if

C     Set the corresponding alabb argument
      if ( ikey .eq. 1 ) then
          checkder = .true.
      else if ( ikey .eq. 2 ) then
          epsgap  = dnum
      else if ( ikey .eq. 3 ) then
          epsfeas = dnum
      else if ( ikey .eq. 4 ) then
          epsopt = dnum
      else if ( ikey .eq. 5 ) then
          maxoutit = inum
      else if ( ikey .eq. 6 ) then
          iprint = inum
      else if ( ikey .eq. 7 ) then
          ncomp = inum
      end if

C     IIERATE
      go to 100

C     END OF LOOP

C     TERMINATIONS

C     CLOSING SPECIFICATION FILE
 200  continue
      close(20)
      go to 500

C     NO SPECIFICATION FILE
 300  continue
      write(*, 9000)
      write(10,9000)
      go to 500

C     ERROR READING THE SPECIFICATION FILE
 400  continue
      write(*, 9010)
      write(10,9010)
      go to 500

C     PRINTING PARAMETERS VALUES
 500  continue

      write(*, 4000)checkder,epsgap,epsfeas,epsopt,maxoutit,iprint,ncomp
      write(10,4000)checkder,epsgap,epsfeas,epsopt,maxoutit,iprint,ncomp

C     NON-EXECUTABLE STATEMENTS
 1000 format(A80)
 2000 format(BN,I20)
 3000 format(BN,F24.0)
 4000 format(/,' ALABB PARAMETERS:',
     +       /,' checkder = ',     19X,L1,
     +       /,' epsgap   = ',8X,1P,D12.4,
     +       /,' epsfeas  = ',8X,1P,D12.4,
     +       /,' epsopt   = ',8X,1P,D12.4,
     +       /,' maxoutit = ',        I20,
     +       /,' iprint   = ',        I20,
     +       /,' ncomp    = ',        I20)
 9000 format(/,' Specification file alabb.dat was not found in the',
     +         ' current directory.',/,' Using default values.')
 9005 format(/,' Specification file alabb.dat is being used.')
 9010 format(/,' Error reading specification file alabb.dat.')
 9020 format(  ' Ignoring unknown keyword ',A32)
 9030 format(  ' Ignoring incomplete keyword ',A32)
 9040 format(1X,A32)
 9041 format(1X,A32,5X,I20)
 9042 format(1X,A32,1X,1P,D24.8)

      end

C     *****************************************************************
C     *****************************************************************
      subroutine checkd(m,n,l,u,wi,wd1,wd2,wd3,wd4,x,macheps,inform)

      implicit none

C     This subrotutine checks the user supplied first and second
C     derivatives subroutines (evalg, evalh, evaljac and evalhc) for
C     computing the objective function gradient and Hessian and the
C     constraints gradients and Hessians, respectively.

C     SCALAR ARGUMENTS
      integer inform,m,n
      double precision macheps

C     ARRAY ARGUMENTS
      integer wi(n)
      double precision l(n),u(n),wd1(n),wd2(n),wd3(n),wd4(n),x(n)

C     Parameters of the subroutine:
C     =============================
C
C     On Entry:
C     =========
C
C     M integer: Number of constraints.
C     ---------
C
C     N integer: Number of variables.
C     ---------
C
C     L double precision l(n): Lower bounds.
C     -----------------------
C
C     U double precision u(n): Upper bounds.
C     -----------------------
C
C     WI double precision wi(n): N-dimensional integer working vector.
C     -----------------------
C
C     WD1 ... WD4, X double precision wd1(n) ... wd4(n), x(n) 
C     -------------------------------------------------------
C
C     N-dimensional double precision working vectors.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     On Return:
C     ==========
C
C     INFORM integer: Output flag.
C     --------------

C     LOCAL SCALARS
      character answer
      integer i,j
      double precision drand,seed,smalll,smallu

      inform = 0

C     SET A RANDOM POINT 

      seed = 17325.0d0
      do i = 1,n
          smalll = max( l(i), - 1.0d0 )
          smallu = min( u(i), smalll + 1.0d0 )
          x(i) = smalll + ( smallu - smalll ) * drand(seed)
      end do
 
      write(*,900)

      do i = 1,n
          write(*,910) i,x(i)
      end do

C     CHECK OBJECTIVE FUNCTION GRADIENT

      write(*,920)

      read(*,*) answer

      if ( answer .eq. 'A' .or. answer .eq. 'a' ) then
          go to 500

      else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then
          call checkg(n,x,wd1,macheps,inform)
          if ( inform .lt. 0 ) then
              return
          end if 
      end if

C     CHECK CONSTRAINTS GRADIENTS

      j = 1

 110  if ( j .le. m ) then

          write(*,930) j

          read(*,*) answer

          if ( answer .eq. 'A' .or. answer .eq. 'a' ) then
              go to 500

          else if ( answer .eq. 'S' .or. answer .eq. 's' ) then
              go to 200

          else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then
              call checkjac(n,x,j,wi,wd1,wd2,macheps,inform)
              if ( inform .lt. 0 ) then
                  return
              end if 
          end if

          j = j + 1

          go to 110

      end if

C     CHECK HESSIAN OF THE OBJECTIVE FUNCTION

 200  continue

      write(*,940)

      read(*,*) answer

      if ( answer .eq. 'A' .or. answer .eq. 'a' ) then
          go to 500

      else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then
          call checkh(n,x,wd1,wd2,wd3,macheps,inform)
          if ( inform .lt. 0 ) then
              return
          end if 
      end if

C     CHECK HESSIANS OF THE CONSTRAINTS

      j = 1

 310  if ( j .le. m ) then

          write(*,950) j

          read(*,*) answer

          if ( answer .eq. 'A' .or. answer .eq. 'a' ) then
              go to 500

          else if ( answer .eq. 'S' .or. answer .eq. 's' ) then
              go to 500

          else if ( answer .eq. 'Y' .or. answer .eq. 'y' ) then
              call checkhc(n,x,j,wi,wd1,wd2,wd3,wd4,macheps,inform)
              if ( inform .lt. 0 ) then
                  return
              end if 
          end if

          j = j + 1

          go to 310

      end if

C     RETURN

 500  continue

      write(*,*) 'Hit any letter to continue.'
      read(*,*) answer

 900  format(/,1X,'Derivatives will be tested at the random point: ')
 910  format(  1X,'x(',I6,') = ',1P,D15.8)
 920  format(/,1X,'Check the gradient of the objective function?',
     +       /,1X,'Type Y(es), N(o) or A(bort checking): ')
 930  format(/,1X,'Check the gradient of constraint ',I5,'?',
     +       /,1X,'Type Y(es), N(o), A(bort checking) or ',
     +            'S(kip constraints gradients): ')
 940  format(/,1X,'Check the Hessian matrix of the objective function?',
     +       /,1X,'Type Y(es), N(o) or A(bort checking): ')
 950  format(/,1X,'Check the Hessian of constraint ',I5,'?',
     +       /,1X,'Type Y(es), N(o), A(bort checking) or ',
     +            'S(kip constraints gradients): ')

      end

C     *****************************************************************
C     *****************************************************************
      subroutine checkg(n,x,g,macheps,inform)

      implicit none

C     This subrotutine checks the user supplied subroutine evalg for 
C     computing the gradient of the objective function using central
C     finite differences with two different discretization steps.

C     SCALAR ARGUMENTS
      integer inform,n
      double precision macheps

C     ARRAY ARGUMENTS
      double precision g(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): Point at which the gradient will be tested.
C     -----------------------
C
C     G double precision g(n) 
C     -----------------------
C
C     N-dimensional double precision working vector.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     On Return:
C     ==========
C
C     INFORM integer: Output flag.
C     --------------

C     LOCAL SCALARS
      integer flag,i
      double precision fminus,fplus,gdiff1,gdiff2,maxerr,step1,step2,tmp

      inform = 0

      call evalg(n,x,g,flag)
      if ( flag .ne. 0 ) then
          inform = - 92
          return
      end if

      write(*,100)

      maxerr = 0.0d0

      do i = 1,n
          tmp  = x(i)

          step1 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d0 )

          x(i) = tmp + step1
          call evalf(n,x,fplus,flag)
          if ( flag .ne. 0 ) then
              inform = - 90
              return
          end if

          x(i) = tmp - step1
          call evalf(n,x,fminus,flag)
          if ( flag .ne. 0 ) then
              inform = - 90
              return
          end if

          gdiff1 = ( fplus - fminus ) / ( 2.0d0 * step1 )

          step2 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d-03 )

          x(i) = tmp + step2
          call evalf(n,x,fplus,flag)
          if ( flag .ne. 0 ) then
              inform = - 90
              return
          end if

          x(i) = tmp - step2
          call evalf(n,x,fminus,flag)
          if ( flag .ne. 0 ) then
              inform = - 90
              return
          end if

          x(i) = tmp

          gdiff2 = ( fplus - fminus ) / ( 2.0d0 * step2 )

          tmp = min( abs( g(i) - gdiff1 ), abs( g(i) - gdiff2 ) )
          write(*,110) i,g(i),gdiff1,gdiff2,tmp
          maxerr = max( maxerr, tmp )

      end do

      write(*,120) maxerr

      return

 100  format(/,1X,'Gradient vector of the objective function.',
     +       /,1X,'Index',13X,'evalg',2X,'Central diff (two different ',
     +            'steps)',4X,'Absolute error')
 110  format(  1X,I5,4(3X,1P,D15.8))
 120  format(  1X,'Maximum absolute error = ',1P,D15.8)

      end

C     *****************************************************************
C     *****************************************************************
      subroutine checkh(n,x,g,gplus1,gplus2,macheps,inform)

      implicit none

C     This subrotutine checks the user supplied subroutine evalh for 
C     computing the Hessian of the objective function using central 
C     finite differences with two different discretization steps.

C     SCALAR ARGUMENTS
      integer inform,n
      double precision macheps

C     ARRAY ARGUMENTS
      double precision g(n),gplus1(n),gplus2(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): Point at which the Hessian will be tested.
C     -----------------------
C
C     G, GPLUS1, GPLUS2 double precision g(n), gplus1(n), gplus2(n) 
C     -------------------------------------------------------------
C
C     N-dimensional double precision working vectors.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     On Return:
C     ==========
C
C     INFORM integer: Output flag.
C     --------------

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      integer flag,i,j,nnzh
      double precision elem,hdiff1,hdiff2,step1,step2,tmp

C     LOCAL ARRAYS
      integer hlin(nmax**2),hcol(nmax**2)
      double precision H(nmax,nmax),hval(nmax**2),maxerr(nmax)

      inform = 0

C     Test viability

      if ( n .gt. nmax ) then
          write(*,*) 'Subroutine CheckH uses dense matrices up to ',
     +               'dimension ',nmax,' times ',nmax,'. The Hessian ',
     +               'checking will be skipped.'
          go to 500
      end if

C     Compute the gradient of the objective function at x

      call evalg(n,x,g,flag)
      if ( flag .ne. 0 ) then
          inform = - 92
          return
      end if

C     Compute the Hessian of the objective function at x and save in a
C     dense matrix

      call evalh(n,x,hlin,hcol,hval,nnzh,flag)
      if ( flag .ne. 0 ) then
          inform = - 94
          return
      end if

      do j = 1,n
          do i = 1,n
              H(i,j) = 0.0d0
          end do
      end do

      do i = 1,nnzh
          H(hlin(i),hcol(i)) = H(hlin(i),hcol(i)) + hval(i)
      end do

C     Test column by column

      write(*,100)

      do j = 1,n

          tmp  = x(j)

          step1 = sqrt( macheps ) * max( abs( tmp ), 1.0d0 )

          x(j) = tmp + step1
          call evalg(n,x,gplus1,flag)
          if ( flag .ne. 0 ) then
              inform = - 92
              return
          end if

          step2 = sqrt( macheps ) * max( abs( tmp ), 1.0d-03 )

          x(j) = tmp + step2
          call evalg(n,x,gplus2,flag)
          if ( flag .ne. 0 ) then
              inform = - 92
              return
          end if

          x(j) = tmp

          write(*,105) j

          maxerr(j) = 0.0d0

          do i = 1,n
              if ( i .ge. j ) then
                  elem = H(i,j)
              else
                  elem = H(j,i)
              end if
              hdiff1 = ( gplus1(i) - g(i) ) / step1
              hdiff2 = ( gplus2(i) - g(i) ) / step2
              tmp = min( abs( elem - hdiff1 ), abs( elem - hdiff2 ) )
              if ( elem .ne. 0.0d0 .or. 
     +        hdiff1 .ne. 0.0d0 .or. hdiff2 .ne. 0.0d0 ) then
                  write(*,110) i,elem,hdiff1,hdiff2,tmp
              end if
              maxerr(j) = max( maxerr(j), tmp )
          end do

          write(*,120) maxerr(j)

      end do

      write(*,*)
      do j = 1,n
          write(*,130) j,maxerr(j)
      end do

 500  continue

      return

 100  format(/,1X,'Hessian matrix of the objective function column by ',
     +            'column.')
 105  format(/,1X,'Column:  ',I6,/,
     +       /,1X,'Index',13X,'evalh',3X,'Incr. Quoc. (two different ',
     +            'steps)',4X,'Absolute error')
 110  format(  1X,I5,4(3X,1P,D15.8))
 120  format(  1X,'Maximum absolute error = ',1P,D15.8)
 130  format(  1X,'Column ',I6,' Maximum absolute error = ',1P,D15.8)

      end

C     *****************************************************************
C     *****************************************************************
      subroutine checkjac(n,x,ind,indjac,valjac,g,macheps,inform)

      implicit none

C     This subrotutine checks the user supplied subroutine evaljac for 
C     computing the gradients of the constraints using central finite 
C     differences with two different discretization steps.

C     SCALAR ARGUMENTS
      integer ind,inform,n
      double precision macheps

C     ARRAY ARGUMENTS
      integer indjac(n)
      double precision g(n),valjac(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) 
C     -----------------------
C
C     Point at which the gradient of the ind-th constraint will be 
C     tested.
C
C     IND integer
C     -----------
C
C     Index of the constraint whose gradient will be tested.
C
C     INDJAC integer indjac(n): N-dimensional integer working vector.
C     ------------------------
C
C     VALJAC, G double precision valjac(n), g(n) 
C     ------------------------------------------
C
C     N-dimensional double precision working vectors.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     On Return:
C     ==========
C
C     INFORM integer: Output flag.
C     --------------

C     LOCAL SCALARS
      integer flag,i,nnzjac
      double precision cminus,cplus,jacdiff1,jacdiff2,maxerr,step1,
     +        step2,tmp

      inform = 0

C     COMPUTE THE GRADIENT OF THE CONSTRAINT AND SAVE IT INTO A DENSE 
C     VECTOR

      call evaljac(n,x,ind,indjac,valjac,nnzjac,flag)
      if ( flag .ne. 0 ) then
          inform = - 93
          return
      end if

      do i = 1,n
          g(i) = 0.0d0
      end do

      do i = 1,nnzjac
          g(indjac(i)) = g(indjac(i)) + valjac(i)
      end do

C     COMPARE WITH CENTRAL FINITE DIFFERENCES

      write(*,100) ind

      maxerr = 0.0d0

      do i = 1,n
          tmp  = x(i)

          step1 = macheps ** (1.0d0/3.0d0) * max( abs( tmp ), 1.0d0 )

          x(i) = tmp + step1
          call evalc(n,x,ind,cplus,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

          x(i) = tmp - step1
          call evalc(n,x,ind,cminus,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

          jacdiff1 = ( cplus - cminus ) / ( 2.0d0 * step1 )

          step2 = macheps * (1.0d0/3.0d0) * max( abs( tmp ), 1.0d-03 )

          x(i) = tmp + step2
          call evalc(n,x,ind,cplus,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

          x(i) = tmp - step2
          call evalc(n,x,ind,cminus,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

          x(i) = tmp

          jacdiff2 = ( cplus - cminus ) / ( 2.0d0 * step2 )

          tmp = min( abs( g(i) - jacdiff1 ), abs( g(i) - jacdiff2 ) )
          if ( g(i) .ne. 0.0d0 .or. 
     +    jacdiff1 .ne. 0.0d0 .or. jacdiff2 .ne. 0.0d0 ) then
              write(*,110) i,g(i),jacdiff1,jacdiff2,tmp
          end if
          maxerr = max( maxerr, tmp )
      end do

      write(*,120) maxerr

      return

 100  format(/,1X,'Gradient vector of constraints ',I5,'.',
     +       /,1X,'Index',11X,'evaljac',2X,'Central diff (two ',
     +            'different steps)',4X,'Absolute error')
 110  format(  1X,I5,4(3X,1P,D15.8))
 120  format(  1X,'Maximum absolute error = ',1P,D15.8)

      end

C     *****************************************************************
C     *****************************************************************
      subroutine checkhc(n,x,ind,indjac,g,gplus1,gplus2,valjac,macheps,
     +inform)

      implicit none

C     This subrotutine checks the user supplied subroutine evalhc for 
C     computing the Hessians of the constraints using finite
C     differences.

C     SCALAR ARGUMENTS
      integer ind,inform,n
      double precision macheps

C     ARRAY ARGUMENTS
      integer indjac(n)
      double precision g(n),gplus1(n),gplus2(n),valjac(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) 
C     -----------------------
C
C     Point at which the Hessian of the ind-th constraint will be 
C     tested.
C
C     IND integer
C     -----------
C
C     Index of the constraint whose gradient will be tested.
C
C     INDJAC integer indjac(n): N-dimensional integer working vector.
C     ------------------------
C
C     G      double precision g(n)
C     GPLUS1 double precision gplus1(n)
C     GPLUS2 double precision gplus2(n)
C     VALJAC double precision valjac(n)
C     ---------------------------------
C
C     N-dimensional double precision working vectors.
C
C     MACHEPS double precision
C     ------------------------
C
C     Smallest positive number such that 1 + macheps is not equal to 1.
C
C     On Return:
C     ==========
C
C     INFORM integer: Output flag.
C     --------------

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      integer flag,i,j,nnzh,nnzjac
      double precision elem,hdiff1,hdiff2,step1,step2,tmp

C     LOCAL ARRAYS
      integer hlin(nmax**2),hcol(nmax**2)
      double precision H(nmax,nmax),hval(nmax**2),maxerr(nmax)

      inform = 0

C     Test viability

      if ( n .gt. nmax ) then
          write(*,*) 'Subroutine CheckHc uses dense matrices up to ',
     +               'dimension ',nmax,' times ',nmax,'. The Hessian ',
     +               'checking will be skipped.'
          go to 500
      end if

C     Compute the constraint gradient at x and save in a dense vector

      call evaljac(n,x,ind,indjac,valjac,nnzjac,flag)
      if ( flag .ne. 0 ) then
          inform = - 93
          return
      end if

      do i = 1,n
          g(i) = 0.0d0
      end do

      do i = 1,nnzjac
          g(indjac(i)) = g(indjac(i)) + valjac(i)
      end do

C     Compute the Hessian of the constraints at x and save in 
C     dense matrix

      call evalhc(n,x,ind,hlin,hcol,hval,nnzh,flag)
      if ( flag .ne. 0 ) then
          inform = - 95
          return
      end if

      do j = 1,n
          do i = 1,n
              H(i,j) = 0.0d0
          end do
      end do

      do i = 1,nnzh
          H(hlin(i),hcol(i)) = H(hlin(i),hcol(i)) + hval(i)
      end do

      write(*,100) ind

      do j = 1,n

          tmp  = x(j)

C         Compute the constraint gradient at xplus1 and save in a 
C         dense vector

          step1 = sqrt( macheps ) * max( abs( tmp ), 1.0d0 )

          x(j) = tmp + step1

          call evaljac(n,x,ind,indjac,valjac,nnzjac,flag)
          if ( flag .ne. 0 ) then
              inform = - 93
              return
          end if

          do i = 1,n
              gplus1(i) = 0.0d0
          end do

          do i = 1,nnzjac
              gplus1(indjac(i)) = valjac(i)
          end do

C         Compute the constraint gradient at xplus2 and save in a 
C         dense vector

          step2 = sqrt( macheps ) * max( abs( tmp ), 1.0d-03 )

          x(j) = tmp + step2

          call evaljac(n,x,ind,indjac,valjac,nnzjac,flag)
          if ( flag .ne. 0 ) then
              inform = - 93
              return
          end if

          do i = 1,n
              gplus2(i) = 0.0d0
          end do

          do i = 1,nnzjac
              gplus2(indjac(i)) = valjac(i)
          end do

          x(j) = tmp

          write(*,105) j

          maxerr(j) = 0.0d0

          do i = 1,n
              if ( i .ge. j ) then
                  elem = H(i,j)
              else
                  elem = H(j,i)
              end if
              hdiff1 = ( gplus1(i) - g(i) ) / step1
              hdiff2 = ( gplus2(i) - g(i) ) / step2
              tmp = min( abs( elem - hdiff1 ), abs( elem - hdiff2 ) )
              if ( elem .ne. 0.0d0 .or.
     +        hdiff1 .ne. 0.0d0 .or. hdiff2 .ne. 0.0d0 ) then
                  write(*,110) i,elem,hdiff1,hdiff2,tmp
              end if
              maxerr(j) = max( maxerr(j), tmp )
          end do

          write(*,120) maxerr(j)

      end do

      write(*,*)
      do j = 1,n
          write(*,130) j,maxerr(j)
      end do

 500  continue

      return

 100  format(/,1X,'Hessian matrix of constraint ',I5,' column by ',
     +            'column.')
 105  format(/,1X,'Column:  ',I6,/,
     +       /,1X,'Index',12X,'evalhc',3X,'Incr. Quoc. (two different ',
     +            'steps)',4X,'Absolute error')
 110  format(  1X,I5,4(3X,1P,D15.8))
 120  format(  1X,'Maximum absolute error = ',1P,D15.8)
 130  format(  1X,'Column ',I6,' Maximum absolute error = ',1P,D15.8)

      end

