!----------------------------------------------------------------------- ! ! 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