mp_pinthread4.f90

Go to the documentation of this file.

!> @file mp_pinthread4.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! This is a CONOPT implementation of the Pindyck model from the
!! GAMS model library. The implementation is similar to the one
!! in pindyck.f90, but this time we first solve the model for
!! a number of periods in parallel, using periods between 21 and 40.
!! We have three loops over periods:
!! First loop is done sequentially
!! The second loop is done in parallel using a static schedule and
!! The third loop is done in parallel using a dynamic schedule.
!! 
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
#if defined(_WIN32) && !defined(_WIN64)
#define dec_directives_win32
#endif
 
Module data
   IMPLICIT NONE
   integer, Parameter :: Tmin = 21, tmax = 60
   integer, Parameter :: Vpp  = 7 ! Variables per period
   integer, Parameter :: epp  = 6 ! Equations per period
   real*8,  Dimension(Tmax, Tmax*Vpp) :: xkeep  ! Store solutions for all periods
   real*8,  Dimension(Tmax, Tmax*Vpp) :: xkeep1 ! Store solutions for all periods
   real*8,  Dimension(Tmax, Tmax*Vpp) :: xkeep2 ! Store solutions for all periods
   Integer, Dimension(Tmax) :: thbase, thstatic, thdynamic, thguided
   Integer :: maxthreadini ! Value of COIGET_MaxThreads at startup
   Logical :: test_fdeval  ! Used to determine if we should test COIDEF_MaxThreads inside FDEval
   Logical :: inparallel   !
   Type controlv
      Integer, Dimension(:), Pointer :: cntvect
   End Type controlv
   Type (controlv), Dimension(:), Pointer :: cntvect1 ! CONOPT Control vectors, one for each thread
!   INTEGER, Dimension(:,:), Allocatable :: CntVect ! Moved to Data to be available in FDEval
End Module data
 
!> Main program. A simple setup and call of CONOPT
!!
Program pinthread
   Use proginfop
   Use conopt
   Use data
   Use omp_lib
   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 :: pin_status     ! Callback for displaying solution status
   Integer, External :: pin_solution   ! Specialized callback for displaying solution values
   Integer, External :: pin_message    ! Specialized Callback for managing messages
   Integer, External :: pin_errmsg     ! Specialized Callback for managing error messages
   Integer, External :: pin_option     ! Define options using a callback. Files cannot be used
                                       ! for multithreading since they can be shared and not
                                       ! read in the proper order.
#ifdef dec_directives_win32
!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 :: Pin_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Option
#endif
!
!  Control vector
!
   Integer, Dimension(:), pointer :: cntvect ! CONOPT Control vector for current thread
   real*8,  Dimension(:), Allocatable :: usrmem
   INTEGER :: coi_error
!
!  Data needed to manage the dynamics and the threads.
!
   Integer :: t, i
   Integer :: maxthreadomp, maxthreadcon, maxthread, thread
   real*8 time0, time1, time2, time3, time4
   Logical dif12, dif13, dif23
   Integer :: ndif12, ndif13, ndif23
 
   call startup
!
!  Create and initialize one Control Vector for each possible thread
!  plus a Usrmem vector to be used to communicate the number of periods
!  to the model routines.
!
   If ( debug ) Then
      open(12,file='Thread1.txt',status='Unknown');
      open(13,file='Thread2.txt',status='Unknown');
      open(14,file='Thread3.txt',status='Unknown');
      open(15,file='Thread4.txt',status='Unknown');
   Endif
!
!  Make a quick test of the COIGET_MaxThread routines. We need a
!  control vector so we create a temporary one.
!
   coi_error = coi_create( cntvect )
   maxthreadomp = omp_get_max_threads()
   write(*,*) 'Initial MaxThread=',maxthreadomp,' from OpenMP'
   maxthreadcon = coiget_maxthreads( cntvect )
   write(*,*) 'Initial MaxThread=',maxthreadcon,' from COIGET_MaxThreads'
   If ( maxthreadcon /= maxthreadomp ) Then
      call flog( "Error in COIGET_MaxThread functions", 1 )
   Endif
   If ( maxthreadomp < 4 ) then
      maxthread = 4
      call omp_set_num_threads(maxthread)
      write(*,*) 'Revised MaxThread=',maxthread
   Else
      maxthread = maxthreadomp
   Endif
   maxthreadcon = coiget_maxthreads( cntvect )
   write(*,*) 'Revised MaxThread=',maxthreadcon,' from COIGET_MaxThreads'
   maxthreadini = maxthreadcon
   coi_error = coi_free( cntvect )
   test_fdeval  = .true.
   inparallel   = .false.
!
!  Create a vector of Control Vectors in CntVect1 and allocate and
!  initialize the individual vectors for each thread.
!
   Allocate( cntvect1( maxthread ) )
   Allocate( usrmem(2*maxthread) )
   DO thread = 1, maxthread
!
!  Creating the control vector for each thread. This will populate the elements
!  in the control vector with information that is independent of the number of periods.
!
      coi_error = max( coi_error, coi_create( cntvect ) )
      cntvect1(thread)%CntVect => cntvect
!
!  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 ) )  ! Objective is constraint 1
!
!  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, pin_status     ) )
      coi_error = max( coi_error, coidef_solution( cntvect, pin_solution   ) )
      coi_error = max( coi_error, coidef_message( cntvect, pin_message    ) )
      coi_error = max( coi_error, coidef_errmsg( cntvect, pin_errmsg     ) )
      If ( debug ) &
      coi_error = max( coi_error, coidef_option( cntvect, pin_option   ) )
!
!  Add the Usrmem address that depends on the thread
!
      coi_error = max( coi_error, coidef_usrmem( cntvect, usrmem(2*thread-1) ) )
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
      coi_error = max( coi_error, coidef_license( cntvect, conopt_license_int_1, conopt_license_int_2, conopt_license_int_3, conopt_license_text) )
#endif
   Enddo
 
   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
!
!  Store special values in Xkeep.
!
   xkeep(tmin:tmax,1:tmax*vpp) = -1234.5d0  ! Store solutions for all periods
!
!  Now perform the actual optimizations, looping over the number of periods
!  In the first loop we run the models sequentially to have a base for
!  comparison.
!
   thread = 1
   cntvect => cntvect1(thread)%CntVect
   time0 = omp_get_wtime()
   DO t = tmin, tmax
      thbase(t) = thread
      write(*,*) 'Sequential: Solving period',t,' using thread',thread
      write(11,*) 'Sequential: Starting period',t,' using thread',thread
      If ( debug ) &
      write(10,*) 'Sequential: Starting period',t,' using thread',thread
      usrmem(2*thread-1) = dfloat(t)      ! Define the content of Usrmem as T
      usrmem(2*thread  ) = dfloat(thread) ! and Thread
!
!  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 ) )
!
      coi_error = coi_solve( cntvect )
      If ( coi_error .ne. 0 ) THEN
         write(*,*)
         write(*,*) '**** Fatal Error while Solving for T=',t,' COI_Error=',coi_error
         write(*,*)
         call flog( "Errors encountered during Solve.", 1 )
      Endif
   Enddo
   time1 = omp_get_wtime() - time0
   xkeep1(tmin:tmax,1:tmax*vpp) = xkeep(tmin:tmax,1:tmax*vpp)
   time0 = omp_get_wtime()
   write(*,*) 'Finished single-thread loop. Going to start static multi-threading loop'
   write(11,*) 'Finished single-thread loop. Going to start static multi-threading loop'
   If ( debug ) &
   write(10,*) 'Finished single-thread loop. Going to start static multi-threading loop'
!
!  Now run the loop again this time in reverse order and in parallel - static.
!
   test_fdeval = .true.
   inparallel  = .true.
   coi_error   = 0
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(Thread, MaxThreadCON, CntVect) SCHEDULE(Static) Reduction(+:COI_Error)
   DO t = tmax, tmin, -1
      thread = 1+omp_get_thread_num()
      thstatic(t) = thread
      cntvect => cntvect1(thread)%CntVect
      write(*,*) 'Static: Solving period',t,' using thread',thread
      write(11,*) 'Static: Starting period',t,' using thread',thread
      If ( debug ) &
      write(10,*) 'Static: Starting period',t,' using thread',thread
      usrmem(2*thread-1) = dfloat(t)      ! Define the content of Usrmem as T
      usrmem(2*thread  ) = dfloat(thread) ! and Thread
!
!  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 ) )
      maxthreadcon = coiget_maxthreads( cntvect )
      If ( maxthreadcon /= 1 ) Then
         coi_error = 1000
         write(*,*)
         write(*,*) '**** Fatal Error in MaxThread while Solving Parallel for T=',t
         write(*,*)
      Else
         coi_error = coi_solve( cntvect )
         coi_error = abs(coi_error)
      Endif
   Enddo
!$OMP END PARALLEL
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while Solving Static OMP loop. COI_Error=',coi_error
      write(*,*)
      call flog( "Errors encountered during Solvein Static OMP loop.", 1 )
   Endif
   time2 = omp_get_wtime() - time0
   xkeep2(tmin:tmax,1:tmax*vpp) = xkeep(tmin:tmax,1:tmax*vpp)
   time0 = omp_get_wtime()
   write(*,*) 'Finished static multi-thread loop. Going to start dynamic multi-threading loop'
   write(11,*) 'Finished static multi-thread loop. Going to start dynamic multi-threading loop'
   If ( debug ) &
   write(10,*) 'Finished static multi-thread loop. Going to start dynamic multi-threading loop'
!
!  Now run the loop again this time in reverse order and in parallel - dynamic.
!
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(Thread, CntVect) SCHEDULE(Dynamic) Reduction(+:COI_Error)
   DO t = tmax, tmin, -1
      thread = 1+omp_get_thread_num()
      thdynamic(t) = thread
      cntvect => cntvect1(thread)%CntVect
      write(*,*) 'Dynamic: Solving period',t,' using thread',thread
      write(11,*) 'Dynamic: Starting period',t,' using thread',thread
      If ( debug ) &
      write(10,*) 'Dynamic: Starting period',t,' using thread',thread
      usrmem(2*thread-1) = dfloat(t)      ! Define the content of Usrmem as T
      usrmem(2*thread  ) = dfloat(thread) ! and Thread
!
!  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 ) )
!
      coi_error = coi_solve( cntvect )
      coi_error = abs(coi_error)
   Enddo
!$OMP END PARALLEL
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while Solving Dynamic OMP loop. COI_Error=',coi_error
      write(*,*)
      call flog( "Errors encountered during Solvein Dynamic OMP loop.", 1 )
   Endif
   time3 = omp_get_wtime() - time0
   time0 = omp_get_wtime()
   write(*,*) 'Finished dynamic multi-thread loop. Going to start guided multi-threading loop'
   write(11,*) 'Finished dynamic multi-thread loop. Going to start guided multi-threading loop'
   If ( debug ) &
   write(10,*) 'Finished dynamic multi-thread loop. Going to start guided multi-threading loop'
!
!  Now run the loop again this time in reverse order and in parallel.
!
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(Thread, CntVect) SCHEDULE(Guided) Reduction(+:COI_Error)
   DO t = tmax, tmin, -1
      thread = 1+omp_get_thread_num()
      thguided(t) = thread
      cntvect => cntvect1(thread)%CntVect
      write(*,*) 'Guided: Solving period',t,' using thread',thread
      write(11,*) 'Guided: Starting period',t,' using thread',thread
      If ( debug ) &
      write(10,*) 'Guided: Starting period',t,' using thread',thread
      usrmem(2*thread-1) = dfloat(t)      ! Define the content of Usrmem as T
      usrmem(2*thread  ) = dfloat(thread) ! and Thread
!
!  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 ) )
!
      coi_error = coi_solve( cntvect )
      coi_error = abs(coi_error)
   Enddo
!$OMP END PARALLEL
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while Solving Guided OMP loop. COI_Error=',coi_error
      write(*,*)
      call flog( "Errors encountered during Solvein Guided OMP loop.", 1 )
   Endif
   write(*,*) 'Finished guided multi-thread loop.'
   write(11,*) 'Finished guided multi-thread loop.'
   If ( debug ) &
   write(10,*) 'Finished guided multi-thread loop.'
   time4 = omp_get_wtime() - time0
!
!  Clean up
!
   DO thread = 1, maxthread
      coi_error = max( coi_error, coi_free( cntvect1(thread)%CntVect ) )
   Enddo
   Deallocate( cntvect1 )
 
   write(*,*)
   write(*,*) 'End of PinThread Model.'
   write(*,*)
!
   dif12 = .false.; ndif12 = 0
   Do t = tmin, tmax
      Do i=1,t*vpp
        if ( xkeep1(t,i) .ne. xkeep2(t,i) ) then
           dif12 = .true.
           ndif12 = ndif12 + 1
           if ( ndif12 <= 10 ) &
              write(*,1000) t, i, thbase(t),'  XBase =',xkeep1(t,i), thstatic(t),'  XStat =',xkeep2(t,i), (xkeep1(t,i)-xkeep2(t,i))
        endif
      Enddo
   Enddo
   write(*,*)
   write(*,*) 'Difference between base and static solutions   =',dif12
   write(*,*)
   dif13 = .false.; ndif13 = 0
   Do t = tmin, tmax
      Do i=1,t*vpp
        if ( xkeep(t,i) .ne. xkeep1(t,i) ) then
           dif13 = .true.
           ndif13 = ndif13 + 1
           if ( ndif13 <= 10 ) &
              write(*,1000) t, i, thbase(t),'  XBase =',xkeep1(t,i), thdynamic(t),'  XDynm =',xkeep(t,i), (xkeep1(t,i)-xkeep(t,i))
        endif
      Enddo
   Enddo
   write(*,*)
   write(*,*) 'Difference between base and dynamic solutions  =',dif13
   write(*,*)
   dif23 = .false.; ndif23 = 0
   Do t = tmin, tmax
      Do i=1,t*vpp
        if ( xkeep(t,i) .ne. xkeep2(t,i) ) then
           dif23 = .true.
           ndif23 = ndif23 + 1
           if ( ndif23 <= 10 ) &
              write(*,1000) t, i, thstatic(t),'  XStat =',xkeep2(t,i), thdynamic(t),'  XDynm =',xkeep(t,i), (xkeep2(t,i)-xkeep(t,i))
        endif
      Enddo
   Enddo
   write(*,*)
   write(*,*) 'Difference between static and dynamic solutions=',dif23
   write(*,*)
1000 Format('T=',i2,' I=',i3,' Thread=',i2,a9,1p,e22.14/ &
                     10x,    ' Thread=',i2,a9,1p,e22.14,'  Diff=',1p,e22.14)
   write(*,*) 'MaxThread =',maxthread
   write(*,*) 'Time         Single thread =',time1
   write(*,*) 'Time Static  Multi thread  =',time2
   write(*,*) 'Time Dynamic Multi thread  =',time3
   write(*,*) 'Time Guided  Multi thread  =',time4
   write(*,*) 'Speedup, static            =',time1/time2
   write(*,*) 'Speedup, dynamic           =',time1/time3
   write(*,*) 'Efficiency, static         =',time1/time2/maxthread
   write(*,*) 'Efficiency, dynamic        =',time1/time3/maxthread
   write(*,*) 'Efficiency, guided         =',time1/time4/maxthread
 
#if defined (nocomp)
!
!  A compiler problem for Intel 11.0.074 and 11.1.065 (in scalar products) means that
!  solutions are not reproducible.
!
   call flog( "Successful Solve", 0 )
#else
   if ( dif12 ) then
      call flog("Base and Static solutions are not identical",1)
   else if ( dif13 ) then
      call flog("Base and Dynamic solutions are not identical",1)
   else if ( dif23 ) then
      call flog("Static and Dynamic solutions are not identical",1)
   else
      call flog( "Successful Solve", 0 )
   endif
#endif
 
end Program pinthread
 
 
!>  Sets runtime options
!!
!! @include{doc} option_params.dox
Integer Function pin_option( ncall, rval, ival, lval, usrmem, name )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Option
#endif
   integer  ncall, ival, lval
   character(Len=*) :: name
   real*8   rval
   real*8   usrmem(*)                                  ! optional user memory
 
   Select Case (ncall)
   case (1)
      name = 'locon3'
      ival = 0
   case (2)
      name = 'lotria'
      ival = 0
   case default
      name = ' '
   end Select
   pin_option = 0
 
end Function pin_option
!
! =====================================================================
!  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 )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_ReadMatrix
#endif
   Use data
      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
   Integer :: t ! Number of periods for this solve
!
!  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.
!
   t = usrmem(1)
   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
!
!  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 )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_FDEval
#endif
   Use data
   Use proginfop
   Use conopt
   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
   Integer :: maxthreads
 
   if ( test_fdeval ) then
!
!  This test is a little artificial and needs a control vector that is not
!  really available here -- but anyway, just test the return value in
!  this case.
!
      maxthreads = coiget_maxthreads( cntvect1(1+thread)%CntVect )
      if ( inparallel ) then
         If ( maxthreads /= 1 ) Then
            write(*,*)
            write(*,*) 'In parallel FDEval. COIGET_MaxThreads=',maxthreads,' and expected 1.'
            write(*,*)
            call flog( "Error in COIGET_MaxThreads in Parallel FDEval", 1 )
         Endif
      Else
         If ( maxthreads /= maxthreadini ) Then
            write(*,*)
            write(*,*) 'In parallel FDEval. COIGET_MaxThreads=',maxthreads,' and expected',maxthreadini
            write(*,*)
            call flog( "Error in COIGET_MaxThreads in Sequential FDEval", 1 )
         Endif
      Endif
      test_fdeval = .false.
   Endif
!
!  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 )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Solution
#endif
 
   Use data
      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(*)
 
   Integer :: t ! number of periods for this solve
 
   t = usrmem(1)
   xkeep(t, 1:n) = xval(1:n) ! Store the primal solution values
   pin_solution = 0
   Return
 
END FUNCTION pin_solution
!
!  Routines for Message, Status, Solution, and ErrMsg.
!
Integer Function pin_status( MODSTA, SOLSTA, ITER, OBJVAL, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Status
#endif
!
   Use proginfop
      IMPLICIT NONE
   INTEGER, Intent(IN)        :: modsta, solsta, iter
   real*8,  Intent(IN)        :: objval
   real*8,  Intent(IN OUT)    :: usrmem(*)
 
   stacalls = stacalls + 1
   obj   = objval
   mstat = modsta
   sstat = solsta
 
   pin_status = 0
 
END Function pin_status
 
Integer Function pin_message( SMSG, DMSG, NMSG, LLEN, USRMEM, MSGV )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_Message
#endif
!
   Use proginfop
      IMPLICIT NONE
 
   INTEGER,            Intent(IN)                 :: smsg, dmsg, nmsg
   INTEGER,            Intent(IN) ,  Dimension(*) :: llen
   CHARACTER(Len=133), Intent(IN),   Dimension(*) :: msgv
   real*8,             Intent(IN OUT)             :: usrmem(*)
 
   integer :: t, thread, i
   t = usrmem(1)
!
! Time-perion Tagged messages in the parallel world
!
   do i = 1, nmsg
      write(11,"(i2,':',A)") t, msgv(i)(1:llen(i))
   enddo
   If ( debug ) then
      thread = usrmem(2)
      do i = 1, dmsg
         write(11+thread,"(i2,':',A)") t, msgv(i)(1:llen(i))
      enddo
   Endif
 
   pin_message = 0
 
END Function pin_message
 
Integer Function pin_errmsg( ROWNO, COLNO, POSNO, MSGLEN, USRMEM, MSG )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Pin_ErrMsg
#endif
!
   Use proginfop
      IMPLICIT NONE
 
   INTEGER,          Intent(IN)                    :: rowno, colno, posno, msglen
   CHARACTER(Len=*), Intent(IN)                    :: msg
   real*8,           Intent(IN OUT)                :: usrmem(*)
!
! Time-period tagged messages in the parallel world
!
   Integer :: t, thread
   t = usrmem(1)
!
!  If Rowno = 0 then the message is about a Column.
!  If Colno = 0 then the message is abuut a Row.
!  Otherwise, the message is about (Rowno,Colno)
!
   IF    ( rowno .EQ. 0 ) THEN  ! Must be a column (error) message
      WRITE(11,1001) t, colno, msg(1:msglen)
   ELSEIF ( colno .EQ. 0 ) THEN  ! Must be a row (error) message
      WRITE(11,1002) t, rowno, msg(1:msglen)
   ELSE                          ! Must be a (row,col) (error) message
      WRITE(11,1000) t, colno, rowno, msg(1:msglen)
   ENDIF
   If ( debug ) Then
      thread = usrmem(2)
      IF    ( rowno .EQ. 0 ) THEN  ! Must be a column (error) message
         WRITE(11+thread,1001) t, colno, msg(1:msglen)
      ELSEIF ( colno .EQ. 0 ) THEN  ! Must be a row (error) message
         WRITE(11+thread,1002) t, rowno, msg(1:msglen)
      ELSE                          ! Must be a (row,col) (error) message
         WRITE(11+thread,1000) t, colno, rowno, msg(1:msglen)
      ENDIF
   Endif
   pin_errmsg = 0
1001  FORMAT(i2,':Variable',i8,' : ',a)
1002  FORMAT(i2,':Equation',i8,' : ',a)
1000  FORMAT(i2,':Variable',i8,' appearing in Equation',i8,' : ',a)
 
END Function pin_errmsg