Log in
!-----------------------------------------------------------------------
!
! Minimise    f
!
! subject to      -0.5 <= x1 <=  INFTY
!                   -2 <= x2 <=  INFTY
!                    0 <= x3 <=  2
!                   -2 <= x4 <=  2
!                         g1 ==  1
!               -INFTY <= g2 <= -1
!                  2.5 <= g3 <=  5
!
! where         f (x1,x2,x3,x4) = x1^2 + 2 x2^2 - x3
!               g1(x1,x2,x3,x4) = x1^2 + x3^2 + x1x3
!               g2(x1,x2,x3,x4) = x3 - x4
!               g3(x1,x2,x3,x4) = x2 + x4
!
! Optimal solution
!                     x*  = (0, 0.5, 1, 2)
!                   f(x*) = -0.5
!                   g(x*) = (1, -1, 2.5)
!
!-----------------------------------------------------------------------
PROGRAM Main
 
! WORHP workspace access macros
#include "worhp_macros.F"
 
  USE Worhp_User ! WORHP user interface module
  USE Worhp_XML  ! module with XML-related functions for WORHP
 
  IMPLICIT NONE
 
  !
  ! WORHP data structures
  !
  ! OptVar contains the variables, constraint values and multipliers.
  ! Workspace encapsulates all workspace and working variables for WORHP.
  ! Params contains all WORHP parameters.
  ! Control contains things for reverse communication flow control.
  !
  TYPE (OptVar)         :: opt
  TYPE (Workspace)      :: wsp
  TYPE (Params)         :: par
  TYPE (Control)        :: cnt
 
  INTEGER(WORHP_INT)    :: status
 
  ! Data structure pointers
  REAL(WORHP_DOUBLE), DIMENSION(:), POINTER :: X, XL, XU, Lambda
  REAL(WORHP_DOUBLE), DIMENSION(:), POINTER :: G, GL, GU, Mu
  INTEGER(WORHP_INT), DIMENSION(:), POINTER :: GType
 
  ! Matrix pointers
  REAL(WORHP_DOUBLE), DIMENSION(:), POINTER :: DFval, DGval, HMval
  INTEGER(WORHP_INT), DIMENSION(:), POINTER :: HMrow, HMcol, DGrow, DGcol, DFrow
 
  ! Check Version of library and header files
  CHECK_WORHP_VERSION
 
  !
  ! Properly zeros everything, or else the following routines
  ! could get confused.
  !
  CALL WorhpPreInit(opt, wsp, par, cnt)
 
  !
  ! Parameter initialisation routine
  ! Must be called before WorhpFromXML, if the optional NoInit flag
  ! is used to initialise all parameters to default values.
  !
  CALL InitParams(status, par)
 
  !
  ! We can now set parameters that may be overruled by those in the
  ! parameter file. This is useful for setting a non-default standard
  ! parameter value that may still be overwritten.
  !
  par%NLPprint = 1   ! Let's prefer the slim output format
                     ! unless the parameter file says differently
 
  !
  ! Parameter XML import routine.
  ! The optional NoInit flag instructs WorhpFromXML not to reset
  ! all parameters to default values (InitParams does this)
  !
  CALL WorhpFromXML(status, "worhp.xml", par = par, NoInit = TRUE)
  IF (status == DataError .OR. status == InitError) STOP
 
  !
  ! WORHP data structure initialisation routine.
  ! Calling this routine prior to WORHP is mandatory.
  ! Before calling WorhpInit, set the problem and matrix dimensions as
  !
  ! opt%N      = number of variables,
  ! opt%M      = number of constraints (lin + nonlin, excluding box con's),
  ! wsp%DF%nnz = nonzero entries of the objective function gradient,
  ! wsp%DG%nnz = nonzero entries of the constraint Jacobian,
  ! wsp%HM%nnz = nonzero entries of the Lagrange Hessian.
  !
  ! Set nnz to 'WorhpMatrix_Init_Dense' to have WorhpInit allocate and
  ! create a dense matrix structure appropriate for the matrix kind and
  ! its dimensions. Setting it to its dense dimension achieves the same.
  !
 
  opt%N = 4  ! This problem has 4 variables
  opt%M = 3  ! and 3 constraints (excluding box constraints)
 
  ! All derivatives for this problem have a sparse structure, so
  ! set the amount of nonzeros here
  wsp%DF%nnz = 3
  wsp%DG%nnz = 6
  wsp%HM%nnz = 1 + opt%N  ! 1 entry on strict lower triangle
                          ! plus full diagonal
 
  CALL WorhpInit(opt, wsp, par, cnt)
  IF (cnt%status /= FirstCall) THEN
     PRINT *, "example: Initialisation failed."
     STOP
  END IF
 
  !
  ! For C interoperability, some arrays have to be C-pointers.
  ! Convert them to Fortran pointers to access them (after WorhpInit)
  !
  CALL C_F_POINTER(opt%X,      X,      [opt%N])
  CALL C_F_POINTER(opt%XL,     XL,     [opt%N])
  CALL C_F_POINTER(opt%XU,     XU,     [opt%N])
  CALL C_F_POINTER(opt%Lambda, Lambda, [opt%N])
 
  CALL C_F_POINTER(opt%G,      G,      [opt%M])
  CALL C_F_POINTER(opt%GL,     GL,     [opt%M])
  CALL C_F_POINTER(opt%GU,     GU,     [opt%M])
  CALL C_F_POINTER(opt%Mu,     Mu,     [opt%M])
  CALL C_F_POINTER(opt%GType,  GType,  [opt%M])
 
  CALL C_F_POINTER(wsp%DF%val, DFval,  [wsp%DF%nnz])
  CALL C_F_POINTER(wsp%DG%val, DGval,  [wsp%DG%nnz])
  CALL C_F_POINTER(wsp%HM%val, HMval,  [wsp%HM%nnz])
 
  !
  ! These pointers give access to the essential user data:
  !
  ! X(1:N)      : Optimisation variables
  ! Lambda(1:N) : Multipliers for the constraints on X ("box constraints")
  ! G(1:M)      : Linear and nonlinear constraints
  ! Mu(1:M)     : Multipliers for the constraints on G
  !
  ! Set initial values of X, Lambda and Mu here.
  ! G need not be initialised.
  !
  X      = [2d0, 2d0, 1d0, 0d0]
  Lambda = ZERO
  Mu     = ZERO
 
  !
  ! Set lower and upper bounds on the variables and constraints.
  ! Use +/-par%Infty to signal "unbounded".
  !
  ! XL and XU are lower and upper bounds ("box constraints") on X.
  ! GL and GU are lower and upper bounds on G.
  !
  XL = [    -5d-1,      -2d0, 0d0, -2d0]
  XU = [par%Infty, par%Infty, 2d0,  2d0]
 
  ! set GL(i) == GU(i) for equality constraints
  GL = [1d0, -par%Infty, 2.5d0]
  GU = [1d0,       -1d0,   5d0]
 
  ! set function types
  opt%FType = WORHP_QUADRATIC
  GType = [WORHP_QUADRATIC, WORHP_LINEAR, WORHP_LINEAR]
 
  !
  ! Define DF as row vector, column index is ommited
  !
  IF (wsp%DF%NeedStructure) THEN
     CALL C_F_POINTER(wsp%DF%row, DFrow,  [wsp%DF%nnz])
     ! only set the nonzero entries, so omit the 4th entry,
     ! which is a structural zero
     DFrow = [1, 2, 3]
  END IF
 
  !
  ! Define DG as CS-matrix
  !
  IF (wsp%DG%NeedStructure) THEN
     CALL C_F_POINTER(wsp%DG%row, DGrow,  [wsp%DG%nnz])
     CALL C_F_POINTER(wsp%DG%col, DGcol,  [wsp%DG%nnz])
     ! only set the nonzero entries in column-major order
     DGrow = [1, 3, 1, 2, 2, 3]
     DGcol = [1, 2, 3, 3, 4, 4]
  END IF
 
  !
  ! Define HM as a diagonal LT-CS-matrix, but only if needed by WORHP
  !
  IF (wsp%HM%NeedStructure) THEN
     CALL C_F_POINTER(wsp%HM%row, HMrow,  [wsp%HM%nnz])
     CALL C_F_POINTER(wsp%HM%col, HMcol,  [wsp%HM%nnz])
     ! only set the nonzero entries, with the strictly lower triangle first,
     ! followed by ALL diagonal entries, so even the entry at (4, 4)
     ! even though it is a structural zero
     HMrow = [3, 1, 2, 3, 4]
     HMcol = [1, 1, 2, 3, 4]
  END IF
 
  !
  ! WORHP Reverse Communication loop.
  ! In every iteration poll GetUserAction for the requested action, i.e. one
  ! of {callWorhp, iterOutput, evalF, evalG, evalDF, evalDG, evalHM, fidif}.
  !
  ! Make sure to reset the requested user action afterwards by calling
  ! DoneUserAction, except for 'callWorhp' and 'fidif'.
  !
  DO WHILE (cnt%status < terminateSuccess .AND. cnt%status > terminateError)
 
     !
     ! WORHP's main routine.
     ! Do not manually reset callWorhp, this is only done by the FD routines.
     !
     IF (GetUserAction(cnt, callWorhp)) THEN
        CALL Worhp(opt, wsp, par, cnt)
        ! Do not reset callWorhp
     END IF
 
     !
     ! Show iteration output.
     ! The call to IterationOutput() may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, iterOutput)) THEN
        CALL IterationOutput(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, iterOutput)
     ENDIF
 
     !
     ! Evaluate the objective function.
     ! The call to UserF may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, evalF)) THEN
        CALL UserF(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, evalF)
     ENDIF
 
     !
     ! Evaluate the constraints.
     ! The call to UserG may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, evalG)) THEN
        CALL UserG(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, evalG)
     ENDIF
 
     !
     ! Evaluate the gradient of the objective function.
     ! The call to UserDF may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, evalDF)) THEN
        CALL UserDF(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, evalDF)
     ENDIF
 
     !
     ! Evaluate the Jacobian of the constraints.
     ! The call to UserDG may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, evalDG)) THEN
        CALL UserDG(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, evalDG)
     ENDIF
 
     !
     ! Evaluate the Hessian matrix of the Lagrange function (L = f + mu*g)
     ! The call to UserHM may be replaced by user-defined code.
     !
     IF (GetUserAction(cnt, evalHM)) THEN
        CALL UserHM(opt, wsp, par, cnt)
        CALL DoneUserAction(cnt, evalHM)
     ENDIF
 
     !
     ! Use finite differences with RC to determine derivatives
     ! Do not reset fidif, this is done by the FD routine.
     !
     IF (GetUserAction(cnt, fidif)) THEN
        CALL WorhpFidif(opt, wsp, par, cnt)
        ! Do not reset fidif
     ENDIF
 
  END DO
 
  !
  ! Translate the WORHP status flag into a meaningful message.
  !
  CALL StatusMsg(opt, wsp, par, cnt)
  !
  ! Deallocate all data structures.
  ! Data structures must not be accessed after this call.
  !
  CALL WorhpFree(opt, wsp, par, cnt)
 
CONTAINS
 
  !
  ! Update F, objective function
  !
  SUBROUTINE UserF(opt, wsp, par, cnt)
    IMPLICIT NONE
    TYPE (OptVar)       :: opt
    TYPE (Workspace)    :: wsp
    TYPE (Params)       :: par
    TYPE (Control)      :: cnt
    INTENT (INOUT)      :: opt, wsp
    INTENT (IN)         :: par, cnt
    opt%F = wsp%ScaleObj * (X(1)**2 + 2d0*X(2)**2 - X(3))
  END SUBROUTINE UserF
 
 
  !
  ! Update G, function of constraints
  !
  SUBROUTINE UserG(opt, wsp, par, cnt)
    IMPLICIT NONE
    TYPE (OptVar)       :: opt
    TYPE (Workspace)    :: wsp
    TYPE (Params)       :: par
    TYPE (Control)      :: cnt
    INTENT (INOUT)      :: opt, wsp
    INTENT (IN)         :: par, cnt
    G = [X(1)**2 + X(3)**2 + X(1)*X(3), X(3) - X(4), X(2) + X(4)]
  END SUBROUTINE UserG
 
 
  !
  ! Update DFval, gradient of objective function
  !
  SUBROUTINE UserDF(opt, wsp, par, cnt)
    IMPLICIT NONE
    TYPE (OptVar)       :: opt
    TYPE (Workspace)    :: wsp
    TYPE (Params)       :: par
    TYPE (Control)      :: cnt
    INTENT (INOUT)      :: opt, wsp
    INTENT (IN)         :: par, cnt
    DFval = wsp%ScaleObj * [2d0*X(1), 4d0*X(2), -ONE]
  END SUBROUTINE UserDF
 
 
  !
  ! Update DGval, Jacobian of constraints
  !
  SUBROUTINE UserDG(opt, wsp, par, cnt)
    IMPLICIT NONE
    TYPE (OptVar)       :: opt
    TYPE (Workspace)    :: wsp
    TYPE (Params)       :: par
    TYPE (Control)      :: cnt
    INTENT (INOUT)      :: opt, wsp
    INTENT (IN)         :: par, cnt
    DGval = [2d0*X(1) + X(3), ONE, X(1) + 2d0*X(3), ONE, -ONE, ONE]
  END SUBROUTINE UserDG
 
 
  !
  ! Update HMval, Hessian of Lagrangian
  !
  SUBROUTINE UserHM(opt, wsp, par, cnt)
    IMPLICIT NONE
    TYPE (OptVar)       :: opt
    TYPE (Workspace)    :: wsp
    TYPE (Params)       :: par
    TYPE (Control)      :: cnt
    INTENT (INOUT)      :: opt, wsp
    INTENT (IN)         :: par, cnt
 
    ! HMval may be reallocated during first WORHP call
    CALL C_F_POINTER(wsp%HM%val, HMval,  [wsp%HM%nnz])
 
    ! Only scale the F part of HM
    HMval = [Mu(1), wsp%ScaleObj * 2d0 + 2d0*Mu(1), wsp%ScaleObj * 4d0, 2d0*Mu(1), ZERO]
 
  END SUBROUTINE UserHM
 
END PROGRAM Main