PROGRAM tminim ! Use minim to maximize the objective function: ! 1/{1 + (x-y)^2} + sin(pi.y.z/2) + exp[-{(x+z)/y - 2}^2] ! with respect to x, y and z. ! We will actually minimize its negative. ! The maximum occurs at x = y = z = +/-sqrt(4n+1) for any integer n. USE Nelder_Mead IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp) :: object, p(3), simp, step(3), stopcr, var(3) INTEGER :: ier, iprint, iquad, maxf, nloop, nop LOGICAL :: first INTERFACE SUBROUTINE objfun(p, func) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp), INTENT(IN) :: p(:) REAL (dp), INTENT(OUT) :: func END SUBROUTINE objfun END INTERFACE ! Set up starting values & step sizes. ! Parameters p(1), p(2), p(3) are x, y & z. p(1) = 0._dp p(2) = 1._dp p(3) = 2._dp step = 0.4_dp nop = 3 ! Set max. no. of function evaluations = 250, print every 10. maxf = 250 iprint = 10 ! Set value for stopping criterion. Stopping occurs when the ! standard deviation of the values of the objective function at ! the points of the current simplex < stopcr. stopcr = 1.d-04 nloop = 6 ! Fit a quadratic surface to be sure a minimum has been found. iquad = 1 ! As function value is being evaluated in REAL (dp), it ! should be accurate to about 15 decimals. If we set simp = 1.d-6, ! we should get about 9 dec. digits accuracy in fitting the surface. simp = 1.d-6 ! Now call MINIM to do the work. first = .true. DO CALL minim(p, step, nop, object, maxf, iprint, stopcr, nloop, & iquad, simp, var, objfun, ier) ! If ier > 0, try a few more function evaluations. IF (ier == 0) EXIT IF (.NOT. first) STOP first = .false. maxf = 100 END DO ! Successful termination. WRITE(*, 900) object, p 900 FORMAT(' Success !'/' Objective function = ', f12.6/ ' at: ', 3f12.6) WRITE(*, 910) var 910 FORMAT(' Elements of var = ', 3f12.6) STOP END PROGRAM tminim SUBROUTINE objfun(p, func) ! This is the subroutine which the user must write. ! Remember that we are minimizing the negative of the function we ! really want to maximize. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp), INTENT(IN) :: p(:) REAL (dp), INTENT(OUT) :: func ! Local variables REAL (dp), PARAMETER :: half = 0.5_dp, one = 1._dp, pi = 3.14159265357989_dp, & two = 2._dp REAL (dp) :: x, y, z x = p(1) y = p(2) z = p(3) func = -one/(one + (x-y)**2) - SIN(half*pi*y*z) - EXP(-((x+z)/y - two)**2) RETURN END SUBROUTINE objfun