fvboth2.f90

Go to the documentation of this file.

!> @file fvboth2.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! This is the pindyck example rewritten to use both the FVforAll
!! and FVincLin way of defining nonlinear functions. All functions, linear as
!! well as nonlinear, will be called, and the functions must return the
!! sum of the linear and nonlinear terms.
!! Similar to FvBoth but with Linear Jacobian elements defined in FDEval
!! 
!!
!! 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
 
!> Main program. A simple setup and call of CONOPT
!!
Program fvboth2
   Use proginfo
   Use conopt
   Use data_t
   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
#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 :: 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
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
!
!  Other variables
!
   INTEGER :: major, minor, patch
 
   call startup
!
!  Create and initialize a Control Vector
!
   coi_error = coi_create( cntvect )
 
!  Write which version of CONOPT we are using.
 
   call coiget_version( major, minor, patch )
   write(*,"('Solving Pindyck Model using CONOPT version ',i2,'.',i2,'.',i2)") major, minor, patch
!
!  Define the number of time periods, T.
!
   t = 16
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
!  Number of variables (excl. slacks): 7 per period
!
   coi_error = max( coi_error, coidef_numvar( cntvect, 7 * t ) )
!
!  Number of equations: 1 objective + 6 per period
!
   coi_error = max( coi_error, coidef_numcon( cntvect, 1 + 6 * 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 ) )
!
!  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 ) )
!
!  Define that all functions should be called.
!
   coi_error = max( coi_error, coidef_fvforall( cntvect, +1 ) )
!
!  Define that functions will include linear terms
!
   coi_error = max( coi_error, coidef_fvinclin( cntvect, +1 ) )
!
!  Turn function debugging on in the initial point to check if it is consistent
!  with FVforAll:
!
   coi_error = max( coi_error, coidef_debugfv( cntvect, -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, 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     ) )
 
#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
 
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping Solve due to setup errors", 1 )
   ENDIF
!
!  Save the solution so we can check the duals:
!
   do_allocate = .true.
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of FVboth2 Model. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-1170.4863d0 ) > 0.0001d0 ) then
      call flog( "Incorrect objective returned", 1 )
   Else
      Call checkdual( 'FVboth2', maximize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
end Program fvboth2
!
! =====================================================================
!  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_t
   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
                                                       ! (not defined here)
   integer, intent (out),    dimension(m)    :: type   ! vector of equation types
   integer, intent (in out), dimension(m)    :: esta   ! vector of initial equation status
                                                       ! (not defined here)
   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
!
!  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 = 7*(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-7) - 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_t
   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
!
   pin_fdeval = 0
   it = (rowno+4) / 6
   is = 7*(it-1)
   if ( rowno == 1 ) then
!
!  This is the linear objective
!  Obj  : sum(t,(1+r)**(1-t)*rev(t)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = 0.0d0
         do it = 1, t
            is = 7*(it-1)
            g = g + x(is+7) * 1.05d0**(1-it)
         enddo
      endif
      if ( mode == 2 .or. mode == 3 ) then
         do it = 1, t
            is = 7*(it-1)
            jac(is+7) = 1.05d0**(1-it)
         enddo
      endif
   else if ( rowno == (it-1)*6+2 ) then
!
!  tddef equation. All terms: td + 0.13*p -0.87*td(t-1)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = x(is+1) + 0.13d0*x(is+6)
         if ( it > 1 ) g = g - 0.87d0*x(is+1-7)
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+1) = 1.0d0
         jac(is+6) = 0.13d0
         if ( it > 1 ) jac(is+1-7) = -0.87d0
      endif
   else if ( rowno == (it-1)*6+3 ) then
!
!  sdef equation. All terms: s - 0.75*s(t-1) -(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 = x(is+3) -h1*h2
         if ( it > 1 ) g = g - 0.75d0*x(is+3-7) ! 0.75*s(t-1) only it t > 1
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+2) = h1*h2*log(1.02d0)/7.d0
         jac(is+6) = -h2*0.1d0
         jac(is+3) = +1.0d0
         if ( it > 1 ) jac(is+3-7) = -0.75d0
      endif
   else if ( rowno == (it-1)*6+4 ) then
!
!  csdef equation. All terms: cs - s - cs(t-1)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = x(is+2) - x(is+3)
         if ( it > 1 ) g = g - x(is+2-7)
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+2) = 1.0d0
         jac(is+3) = -1.0d0
         if ( it > 1 ) jac(is+2-7) = -1.0d0
      endif
   else if ( rowno == (it-1)*6+5 ) then
!
!  ddef equation. All terms: -td + s + d
!
      if ( mode == 1 .or. mode == 3 ) then
         g = -x(is+1) + x(is+3) + x(is+4)
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+1) = -1.0d0
         jac(is+3) =  1.0d0
         jac(is+4) =  1.0d0
      endif
   else if ( rowno == (it-1)*6+6 ) then
!
!  rdef equation. All terms: d + r - r(t-1)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = x(is+4) + x(is+5)
         if ( it > 1 ) g = g - x(is+5-7)
      endif
      if ( mode == 2 .or. mode == 3 ) then
         jac(is+4) = 1.0d0
         jac(is+5) = 1.0d0
         if ( it > 1 ) jac(is+5-7) = -1.0d0
      endif
   elseif ( rowno == (it-1)*6+7 ) then
!
!  revdef equation. All terms: rev - d*(p-250/r)
!
      if ( mode == 1 .or. mode == 3 ) then
         g = x(is+7) -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)
         jac(is+7) = +1.0d0
      endif
   else
!
!  We should have covered everything above so this is an error.
!
      pin_fdeval = 1
   endif
 
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
!
!  Specialized solution callback routine with names for variables and constraints
!
   Use proginfo
   Use data_t
   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' /)
 
   WRITE(10,"(/' Variable   Solution value    Reduced cost    B-stat'/)")
   i = 0
   do it = 1, t
      DO i1 = 1, 7
         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, 6
         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
 
   solcalls = solcalls + 1
   pin_solution = 0
 
END Function pin_solution