pinadd2.f90

Go to the documentation of this file.

!> @file pinadd2.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! This is a CONOPT implementation of the Pindyck model from the
!! GAMS model library. The implementation is similar to the one
!! in pindyck, but this time we first solve the model for 16
!! periods and then we gradually increase the number of periods
!! one at a time up to 20.
!! 
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
Module pinadd2data
   integer :: t
   integer, Parameter :: tmin = 16, tmax = 20
   integer, Parameter :: vpp  = 7 ! Variables per period
   integer, Parameter :: epp  = 6 ! Equations per period
   Integer, Parameter :: hpp  = 5 ! Hessian elements per period -- see in 2DLagr* routines
   real*8,  Dimension(Tmax*Vpp) :: xkeep  ! Store previous solutions
   Integer, Dimension(Tmax*Vpp) :: xstas  ! Store previous status
   Integer, Dimension(Tmax*Epp) :: estat  ! for variables and equations
End Module pinadd2data
 
!> Main program. A simple setup and call of CONOPT
!!
Program pinadd2
   Use proginfo
   Use coidef
   Use pinadd2data
   Implicit none
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: pin_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: pin_fdeval     ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: pin_solution   ! Specialized callback for displaying solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
   Integer, External :: pin_2dlagrstr  ! Hessian structure evaluation routine
   Integer, External :: pin_2dlagrval  ! Hessian value evaluation routine
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_2DLagrVal
#endif
!
!  Control vector
!
   INTEGER :: numcallback
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   call startup
!
!  Create and initialize a Control Vector
!
   numcallback = coidef_size()
   Allocate( cntvect(numcallback) )
   coi_error = coidef_inifort( cntvect )
!
!  Define the number of time periods, T.
!
   t = tmin
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
!  Number of variables (excl. slacks): Vpp per period
!
   coi_error = max( coi_error, coidef_numvar( cntvect, vpp * t ) )
!
!  Number of equations: 1 objective + Epp per period
!
   coi_error = max( coi_error, coidef_numcon( cntvect, 1 + epp * t ) )
!
!  Number of nonzeros in the Jacobian. See the counting in ReadMatrix below:
!  For each period there is 1 in the objective, 16 for unlagged
!  variables and 4 for lagged variables.
!
   coi_error = max( coi_error, coidef_numnz( cntvect, 17 * t + 4 * (t-1) ) )
!
!  Number of nonlinear nonzeros. 5 unlagged for each period.
!
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 5 * t ) )
!
!  Number of nonzeros in the Hessian. 5 for each period.
!
   coi_error = max( coi_error, coidef_numhess( cntvect, hpp * t ) )
!
!  Direction: +1 = maximization.
!
   coi_error = max( coi_error, coidef_optdir( cntvect, 1 ) )
!
!  Objective: Constraint no 1
!
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )
!
!  Turn Hessian debugging on in the initial point
!
   coi_error = max( coi_error, coidef_debug2d( cntvect, -1 ) )
!
!  Options file = test.opt
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'test.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, pin_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, pin_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, pin_solution   ) )
   coi_error = max( coi_error, coidef_message( cntvect, std_message    ) )
   coi_error = max( coi_error, coidef_errmsg( cntvect, std_errmsg     ) )
   coi_error = max( coi_error, coidef_2dlagrstr( cntvect, pin_2dlagrstr  ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, pin_2dlagrval  ) )
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, license_int_1, license_int_2, license_int_3, license_text) )
#endif
 
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping First Solve due to setup errors", 1 )
   ENDIF
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while Solving first model.'
      write(*,*)
      call flog( "Errors encountered during first Solve.", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called during first Solve", 1 )
   else if ( mstat /= 2 .or. sstat /= 1 ) then
      call flog( "Incorrect Model or Solver Status during first solve", 1 )
   ENDIF
!
!  Now increase the time horizon gradually up to 20.
!  All sizes that depend on T must be registered again.
!
   Do t = tmin+1, tmax
      coi_error = max( coi_error, coidef_numvar( cntvect, 7 * t ) )
      coi_error = max( coi_error, coidef_numcon( cntvect, 1 + 6 * t ) )
      coi_error = max( coi_error, coidef_numnz( cntvect, 17 * t + 4 * (t-1) ) )
      coi_error = max( coi_error, coidef_numnlnz( cntvect, 5 * t ) )
      coi_error = max( coi_error, coidef_numhess( cntvect, hpp * t ) )
!
!  This time we have status information copied from a previous solve, i.e.
!  IniStat = 2
!
      coi_error = max( coi_error, coidef_inistat( cntvect, 2 ) )
!
!  Turn Hessian debugging off again
!
      coi_error = max( coi_error, coidef_debug2d( cntvect, 0 ) )
      If ( coi_error .ne. 0 ) THEN
         write(*,*)
         write(*,*) '**** Fatal Error while revising CONOPT Sizes.'
         write(*,*)
         call flog( "Skipping Later Solve due to setup errors", 1 )
      ENDIF
      coi_error = coi_solve( cntvect )
      If ( coi_error .ne. 0 ) THEN
         write(*,*)
         write(*,*) '**** Fatal Error while Solving for T=',t
         write(*,*)
         call flog( "Errors encountered during later Solve.", 1 )
      elseif ( stacalls == 0 .or. solcalls == 0 ) then
         call flog( "Status or Solution routine was not called during later Solve", 1 )
      else if ( mstat /= 2 .or. sstat /= 1 ) then
         call flog( "Incorrect Model or Solver Status during later solve", 1 )
      ENDIF
   Enddo
   write(*,*)
   write(*,*) 'End of Pinadd2 Model. Return code=',coi_error
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
end Program pinadd2
!
! =====================================================================
!  Define information about the model structure
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function pin_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, n, m, nz, usrmem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_ReadMatrix
#endif
   Use pinadd2data
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: m      ! number of constraints
   integer, intent (in)                      :: nz     ! number of nonzeros
   real*8,  intent (in out), dimension(n)    :: lower  ! vector of lower bounds
   real*8,  intent (in out), dimension(n)    :: curr   ! vector of initial values
   real*8,  intent (in out), dimension(n)    :: upper  ! vector of upper bounds
   integer, intent (in out), dimension(n)    :: vsta   ! vector of initial variable status
                                                       ! Defined during restarts
   integer, intent (out),    dimension(m)    :: type   ! vector of equation types
   integer, intent (in out), dimension(m)    :: esta   ! vector of initial equation status
                                                       ! Defined during restarts
   real*8,  intent (in out), dimension(m)    :: rhs    ! vector of right hand sides
   integer, intent (in out), dimension(n+1)  :: colsta ! vector with start of column indices
   integer, intent (out),    dimension(nz)   :: rowno  ! vector of row numbers
   integer, intent (in out), dimension(nz)   :: nlflag ! vector of nonlinearity flags
   real*8,  intent (in out), dimension(nz)   :: value  ! vector of matrix values
   real*8   usrmem(*)                                  ! optional user memory
 
   Integer :: it, is, i, icol, iz, iold
!
!  Define the information for the columns.
!
!  We should not supply status information, vsta.
!
!  We order the variables as follows:
!  td, cs, s, d, r, p, and rev. All variables for period 1 appears
!  first followed by all variables for period 2 etc.
!
!  td, cs, s, and d have lower bounds of 0, r and p have lower
!  bounds of 1, and rev has no lower bound.
!  All have infinite upper bounds (default).
!  The initial value of td is 18, s is 7, cs is 7*t, d is td-s,
!  p is 14, and r is r(t-1)-d. No initial value for rev.
!
   do it = 1, t
      is = vpp*(it-1)
      lower(is+1) = 0.d0
      lower(is+2) = 0.d0
      lower(is+3) = 0.d0
      lower(is+4) = 0.d0
      lower(is+5) = 1.d0
      lower(is+6) = 1.d0
      curr(is+1)  = 18.d0
      curr(is+2)  = 7.d0*it
      curr(is+3)  = 7.d0
      curr(is+4)  = curr(is+1) - curr(is+3)
      if ( it .gt. 1 ) then
         curr(is+5) = curr(is+5-vpp) - curr(is+4)
      else
         curr(is+5) = 500.d0 - curr(is+4)
      endif
      curr(is+6)  = 14.d0
   enddo
   If ( t > tmin ) Then
!
!  This is a restart: Use the initial values from last solve for
!  the variables in the first periods and extrapolate the last
!  period using a linear extrapolation between the last two
!
      iold = vpp*(t-1)
      curr(1:iold) = xkeep(1:iold)
      do i = 1, vpp
         curr(iold+i) = curr(iold+i-vpp) + (curr(iold+i-vpp)-curr(iold+i-2*vpp))
      enddo
!
!  The variables from the last solve are given the old status and
!  the variables in the new period are given those in the last
!  pariod. Similarly with the Equation status:
!
      vsta(1:iold) = xstas(1:iold)
      do i = 1, vpp
         vsta(iold+i) = xstas(iold+i-vpp)
      enddo
      iold = 6*(t-1)+1
      esta(1:iold) = estat(1:iold)
      do i = 1, epp
         esta(iold+i) = estat(iold+i-epp)
      enddo
   endif
!
!  Define the information for the rows
!
!  We order the constraints as follows: The objective is first,
!  followed by tddef, sdef, csdef, ddef, rdef, and revdef for
!  the first period, the same for the second period, etc.
!
!  The objective is a nonbinding constraint:
!
   type(1) = 3
!
!  All others are equalities:
!
   do i = 2, m
      type(i) = 0
   enddo
!
!  Right hand sides: In all periods except the first, only tddef
!  has a nonzero right hand side of 1+2.3*1.015**(t-1).
!  In the initial period there are contributions from lagged
!  variables in the constraints that have lagged variables.
!
   do it = 1, t
      is = 1 + 6*(it-1)
      rhs(is+1) = 1.d0+2.3d0*1.015d0**(it-1)
   enddo
!
!  tddef: + 0.87*td(0)
!
   rhs(2) = rhs(2) + 0.87d0*18.d0
!
!  sdef: +0.75*s(0)
!
   rhs(3) =  0.75d0*6.5d0
!
!  csdef: +1*cs(0)
!
   rhs(4) =  0.d0
!
!  rdef: +1*r(0)
!
   rhs(6) = 500.d0
!
!  Define the structure and content of the Jacobian:
!  To help define the Jacobian pattern and values it can be useful to
!  make a picture of the Jacobian. We describe the variables for one
!  period and the constraints they are part of:
!
!              td   cs   s    d    r    p    rev
!  Obj                                      (1+r)**(1-t)
! Period t:
!  tddef       1.0                      0.13
!  sdef             NL   1.0            NL
!  csdef            1.0 -1.0
!  ddef       -1.0       1.0  1.0
!  rdef                       1.0  1.0
!  revdef                     NL   NL   NL   1.0
! Period t+1:
!  tddef      -0.87
!  sdef                 -0.75
!  csdef           -1.0
!  ddef
!  rdef                           -1.0
!  revdef
!
!  The Jacobian has to be sorted column-wise so we will just define
!  the elements column by column according to the table above:
!
   iz   = 1
   icol = 1
   do it = 1, t
!
!  is points to the position before the first equation for the period
!
      is = 1 + 6*(it-1)
!
!  Column td:
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+1
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+4
      value(iz)  = -1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      if ( it .lt. t ) then
         rowno(iz)  = is+7
         value(iz)  = -0.87d0
         nlflag(iz) = 0
         iz         = iz + 1
      endif
!
!  Column cs
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+2
      nlflag(iz) = 1
      iz         = iz + 1
      rowno(iz)  = is+3
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      if ( it .lt. t ) then
         rowno(iz)  = is+9
         value(iz)  = -1.d0
         nlflag(iz) = 0
         iz         = iz + 1
      endif
!
!  Column s
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+2
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+3
      value(iz)  = -1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+4
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      if ( it .lt. t ) then
         rowno(iz)  = is+8
         value(iz)  = -0.75d0
         nlflag(iz) = 0
         iz         = iz + 1
      endif
!
!  Column d:
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+4
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+5
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+6
      nlflag(iz) = 1
      iz         = iz + 1
!
!  Column r:
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+5
      value(iz)  = +1.d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+6
      nlflag(iz) = 1
      iz         = iz + 1
      if ( it .lt. t ) then
         rowno(iz)  = is+11
         value(iz)  = -1.d0
         nlflag(iz) = 0
         iz         = iz + 1
      endif
!
!  Column p:
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = is+1
      value(iz)  = +0.13d0
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+2
      nlflag(iz) = 1
      iz         = iz + 1
      rowno(iz)  = is+6
      nlflag(iz) = 1
      iz         = iz + 1
!
!  Column rev:
!
      colsta(icol) = iz
      icol       = icol + 1
      rowno(iz)  = +1
      value(iz)  = 1.05d0**(1-it)
      nlflag(iz) = 0
      iz         = iz + 1
      rowno(iz)  = is+6
      value(iz)  = 1.d0
      nlflag(iz) = 0
      iz         = iz + 1
   enddo
   colsta(icol) = iz
 
   pin_readmatrix = 0
 
end Function pin_readmatrix
!
! =====================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function pin_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nz, thread, usrmem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_FDEval
#endif
   Use pinadd2data
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: rowno  ! number of the row to be evaluated
   integer, intent (in)                      :: nz     ! number of nonzeros in this row
   real*8,  intent (in),     dimension(n)    :: x      ! vector of current solution values
   real*8,  intent (in out)                  :: g      ! constraint value
   real*8,  intent (in out), dimension(n)    :: jac    ! vector of derivatives for current constraint
   integer, intent (in),     dimension(nz)   :: jcnm   ! list of variables that appear nonlinearly
                                                       ! in this row. Ffor information only.
   integer, intent (in)                      :: mode   ! evaluation mode: 1 = function value
                                                       ! 2 = derivatives, 3 = both
   integer, intent (in)                      :: ignerr ! if 1 then errors can be ignored as long
                                                       ! as errcnt is incremented
   integer, intent (in out)                  :: errcnt ! error counter to be incremented in case
                                                       ! of function evaluation errors.
   integer, intent (in)                      :: thread
   real*8   usrmem(*)                                  ! optional user memory
 
   integer :: it, is
   real*8  :: h1, h2
!
!  Compute the number of the period
!
   it = (rowno+4) / epp
   is = vpp*(it-1)
   if ( rowno == (it-1)*epp+3 ) then
!
!  sdef equation. Nonlinear term = -(1.1+0.1*p)*1.02**(-cs/7)
!
      h1 = (1.1d0+0.1d0*x(is+6))
      h2 = 1.02d0**(-x(is+2)/7.d0)
      if ( mode == 1 .or. mode == 3 ) then
         g = -h1*h2
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+2) = h1*h2*log(1.02d0)/7.d0
         jac(is+6) = -h2*0.1d0
      endif
   elseif ( rowno == (it-1)*epp+7 ) then
!
!  revdef equation. Nonlinear term = -d*(p-250/r)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = -x(is+4)*(x(is+6)-250.d0/x(is+5))
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+4) = -(x(is+6)-250.d0/x(is+5))
         jac(is+5) = -x(is+4)*250d0/x(is+5)**2
         jac(is+6) = -x(is+4)
      endif
   else
!
!  Error - this equation is not nonlinear
!
      write(*,*)
      write(*,*) 'Error. FDEval called with rowno=',rowno
      write(*,*)
      pin_fdeval = 1
      Return
   endif
   pin_fdeval = 0
 
end Function pin_fdeval
 
Integer Function pin_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Solution
#endif
 
   Use pinadd2data
   Use proginfo
   IMPLICIT NONE
   INTEGER, Intent(IN)                  :: n, m
   INTEGER, Intent(IN),    Dimension(N) :: xbas, xsta
   INTEGER, Intent(IN),    Dimension(M) :: ybas, ysta
   real*8,  Intent(IN),    Dimension(N) :: xval, xmar
   real*8,  Intent(IN),    Dimension(M) :: yval, ymar
   real*8,  Intent(IN OUT)              :: usrmem(*)
   character*6, parameter, dimension(7) :: vname = (/'td    ','cs    ','s     ','d     ','r     ','p     ','rev   '/)
   character*6, parameter, dimension(6) :: ename = (/'tddef ','sdef  ','csdef ','ddef  ','rdef  ','revdef'/)
 
   INTEGER :: i, it, i1
   CHARACTER*5, Parameter, Dimension(4) :: stat = (/ 'Lower','Upper','Basic','Super' /)
 
   solcalls = solcalls + 1
   If ( t < tmax ) Then
      Write(*,"(A,i3)") 'Saving primal solution and status information for T =',t
      xkeep(1:n) = xval(1:n)
      xstas(1:n) = xbas(1:n)
      estat(1:m) = ybas(1:m)
      pin_solution = 0
      Return
   Endif
 
   WRITE(10,"(/' Variable   Solution value    Reduced cost    B-stat'/)")
   i = 0
   do it = 1, t
      DO i1 = 1, vpp
         i = i + 1
         WRITE(10,"(1X,A6,i2,1p,E20.6,E16.6,4X,A5 )") vname(i1), it, xval(i), xmar(i), stat(1+xbas(i))
      ENDDO
   enddo
 
   WRITE(10,"(/' Constrnt   Activity level    Marginal cost   B-stat'/)")
   i = 1
   WRITE(10,"(1x,'Objective',1P,E19.6,E16.6,4X,A5 )") yval(i), ymar(i), stat(1+ybas(i))
   do it = 1, t
      do i1 = 1, epp
         i = i + 1
         WRITE(10,"(1x,A6,i2,1P,E20.6,E16.6,4X,A5 )") ename(i1),it, yval(i), ymar(i), stat(1+ybas(i))
      enddo
   ENDDO
 
   pin_solution = 0
 
END Function pin_solution
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function pin_2dlagrstr( HSRW, HSCL, NODRV, N, M, NHESS, UsrMem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_2DLagrStr
#endif
   USE pinadd2data
   Implicit None
 
   Integer, Intent (IN)                    :: n, m, nhess
   Integer, Intent (IN OUT)                :: nodrv
   Integer, Dimension(NHess), Intent (Out) :: hsrw, hscl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   Integer :: it, is, ic
!
!  There are two nonlinear constraints, Sdef and RevDef and they have
!  the following Hessians (for each time period):
!
!  Revdef equation. Nonlinear term = -d*(p-250/r)
!  Hessian (diagonal and lower triangle): A total of 3 nonzero elements per period
!                4           5           6
!                d           r           p
!  4 : d
!  5 : r     -250/r**2   500*d/r**3
!  6 : p        -1
!
!  Sdef: Nonlinear term = -(1.1+0.1p)*1.02**(-cs/7)
!  Hessian (diagonal and lower triangle): A total of 2 nonzero elements per period
!                           2                                  6
!                          cs                                  p
!  2 : cs  -(1.1+0.1p)*1.02**(-cs/7)*(log(1.02)/7)**2
!  6 : p           0.1*1.02**(-cs/7)*(log(1.02)72)
!
!  The two Hessian matrices do not overlap so the total number of nonzeros is
!  5 per period. The structure using numerical indices for rows and columns and
!  running indices for the Hessian element is:
!
!       2   4   5   6
!  2    1
!  4
!  5        3   5
!  6    2   4
!
!  and the same with names of variables and contributing equations
!       cs      d       r     p
!  cs  sdef
!  d
!  r         revdef  revdef
!  p   sdef  revdef
!
!
!  Define structure of Hessian
!
   if ( n /= vpp*t ) then
      write(*,*) 'Struct: Number of variables should be=',vpp*t,' but CONOPT uses',n
      write(*,*) 'Number of periods=',t
      pin_2dlagrstr = 1; Return
      stop
   endif
   if ( m /= epp*t+1 ) then
      write(*,*) 'Struct: Number of equations should be=',epp*t+1,' but CONOPT uses',m
      write(*,*) 'Number of periods=',t
      pin_2dlagrstr = 1; Return
   endif
   Do it = 1, t
      is = hpp*(it-1)
      ic = vpp*(it-1)
      hscl(is+1) = ic+2; hsrw(is+1) = ic+2
      hscl(is+2) = ic+2; hsrw(is+2) = ic+6
      hscl(is+3) = ic+4; hsrw(is+3) = ic+5
      hscl(is+4) = ic+4; hsrw(is+4) = ic+6
      hscl(is+5) = ic+5; hsrw(is+5) = ic+5
   Enddo
 
   pin_2dlagrstr = 0
 
End Function pin_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function pin_2dlagrval( X, U, HSRW, HSCL, HSVL, NODRV, N, M, NHESS, UsrMem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_2DLagrVal
#endif
   USE pinadd2data
   Implicit None
 
   Integer, Intent (IN)                    :: n, m, nhess
   Integer, Intent (IN OUT)                :: nodrv
   real*8,  Dimension(N), Intent (IN)      :: x
   real*8,  Dimension(M), Intent (IN)      :: u
   Integer, Dimension(Nhess), Intent (IN)  :: hsrw, hscl
   real*8,  Dimension(NHess), Intent (Out) :: hsvl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   Integer :: it, is, ic, ie
   real*8  :: cs, d, r, p
   real*8  :: h1, h2
!
!  There are two nonlinear constraints, Sdef and RevDef and they have
!  the following Hessians (for each time period):
!
!  Revdef equation. Nonlinear term = -d*(p-250/r)
!  Hessian (diagonal and lower triangle): A total of 3 nonzero elements per period
!                4           5           6
!                d           r           p
!  4 : d
!  5 : r     -250/r**2   500*d/r**3
!  6 : p        -1
!
!  Sdef: Nonlinear term = -(1.1+0.1p)*1.02**(-cs/7)
!  Hessian (diagonal and lower triangle): A total of 2 nonzero elements per period
!                           2                                  6
!                          cs                                  p
!  2 : cs  -(1.1+0.1p)*1.02**(-cs/7)*(log(1.02)/7)**2
!  6 : p           0.1*1.02**(-cs/7)*(log(1.02)72)
!
!  The two Hessian matrices do not overlap so the total number of nonzeros is
!  5 per period. The structure using numerical indices for rows and columns and
!  running indices for the Hessian element is:
!
!       2   4   5   6
!  2    1
!  4
!  5        3   5
!  6    2   4
!
!  and the same with names of variables and contributing equations
!       cs      d       r     p
!  cs  sdef
!  d
!  r         revdef  revdef
!  p   sdef  revdef
!
!
!  Normal Evaluation mode
!
   if ( n /= vpp*t ) then
      write(*,*) 'Eval: Number of variables should be=',vpp*t,' but CONOPT uses',n
      write(*,*) 'Number of periods=',t
      pin_2dlagrval = 1; REturn
   endif
   if ( m /= epp*t+1 ) then
      write(*,*) 'Eval: Number of equations should be=',epp*t+1,' but CONOPT uses',m
      write(*,*) 'Number of periods=',t
      pin_2dlagrval = 1; Return
   endif
   Do it = 1, t
      is = hpp*(it-1)
      ic = vpp*(it-1)
      ie = epp*(it-1) + 1
      cs = x(ic+2)
      d  = x(ic+4)
      r  = x(ic+5)
      p  = x(ic+6)
      h1 = 1.1d0+0.1d0*p
      h2 = 1.02d0**(-cs/7.d0)
      hsvl(is+1) = -h1*h2*(log(1.02)/7.d0)**2    * u(ie+2)  ! Sdef = equation 2
      hsvl(is+2) =  0.1d0*h2*(log(1.02)/7.d0)    * u(ie+2)
      hsvl(is+3) = -250.d0/r**2                  * u(ie+6)  ! Revdef = equation 6
      hsvl(is+4) = -1.d0                         * u(ie+6)
      hsvl(is+5) =  500.d0*d/ r**3               * u(ie+6)
   enddo
 
   pin_2dlagrval = 0
 
End Function pin_2dlagrval