LCOV - code coverage report
Current view: top level - src - mo_opt_functions.f90 (source / functions) Hit Total Coverage
Test: forces coverage Lines: 40 971 4.1 %
Date: 2024-03-13 19:03:28 Functions: 3 95 3.2 %

          Line data    Source code
       1             : !> \file mo_opt_functions.f90
       2             : !> \brief \copybrief mo_opt_functions
       3             : !> \details \copydetails mo_opt_functions
       4             : 
       5             : !> \brief Added for testing purposes of test_mo_sce, test_mo_dds, test_mo_mcmc
       6             : !> \author Matthias Cuntz
       7             : !> \date Jul 2012
       8             : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
       9             : !! FORCES is released under the LGPLv3+ license \license_note
      10             : MODULE mo_opt_functions
      11             : 
      12             :   use mo_kind, only: i4, dp
      13             :   use mo_optimization_utils, only: eval_interface
      14             : 
      15             :   IMPLICIT NONE
      16             : 
      17             :   PRIVATE
      18             : 
      19             :   ! ------------------------------------------------------------------
      20             :   ! test_min package of John Burkardt
      21             :   PUBLIC :: quadratic                         ! Simple quadratic, (x-2)^2+1.
      22             :   PUBLIC :: quadratic_exponential             ! Quadratic plus exponential, x^2 + e^(-x).
      23             :   PUBLIC :: quartic                           ! Quartic, x^4 + 2x^2 + x + 3.
      24             :   PUBLIC :: steep_valley                      ! Steep valley, e^x + 1/(100x).
      25             :   PUBLIC :: steep_valley2                     ! Steep valley, e^x - 2x + 1/(100x) - 1/(1000000x^2)
      26             :   PUBLIC :: dying_snake                       ! The dying snake, ( x + sin(x) ) * e^(-x^2).
      27             :   PUBLIC :: thin_pole                         ! The "Thin Pole", x^2+1+log((pi-x)^2)/pi^4
      28             :   PUBLIC :: oscillatory_parabola              ! The oscillatory parabola
      29             :   PUBLIC :: cosine_combo                      ! The cosine combo
      30             :   PUBLIC :: abs1                              ! 1 + |3x-1|
      31             :   ! ------------------------------------------------------------------
      32             :   ! test_opt package of John Burkardt
      33             :   PUBLIC :: fletcher_powell_helical_valley    ! The Fletcher-Powell helical valley function, N = 3.
      34             :   PUBLIC :: biggs_exp6                        ! The Biggs EXP6 function, N = 6.
      35             :   PUBLIC :: gaussian                          ! The Gaussian function, N = 3.
      36             :   PUBLIC :: powell_badly_scaled               ! The Powell badly scaled function, N = 2.
      37             :   PUBLIC :: box_3dimensional                  ! The Box 3-dimensional function, N = 3.
      38             :   PUBLIC :: variably_dimensioned              ! The variably dimensioned function, 1 <= N.
      39             :   PUBLIC :: watson                            ! The Watson function, 2 <= N.
      40             :   PUBLIC :: penalty1                          ! The penalty function #1, 1 <= N.
      41             :   PUBLIC :: penalty2                          ! The penalty function #2, 1 <= N.
      42             :   PUBLIC :: brown_badly_scaled                ! The Brown badly scaled function, N = 2.
      43             :   PUBLIC :: brown_dennis                      ! The Brown and Dennis function, N = 4.
      44             :   PUBLIC :: gulf_rd                           ! The Gulf R&D function, N = 3.
      45             :   PUBLIC :: trigonometric                     ! The trigonometric function, 1 <= N.
      46             :   PUBLIC :: ext_rosenbrock_parabolic_valley   ! The extended Rosenbrock parabolic valley function, 1 <= N.
      47             :   PUBLIC :: ext_powell_singular_quartic       ! The extended Powell singular quartic function, 4 <= N.
      48             :   PUBLIC :: beale                             ! The Beale function, N = 2.
      49             :   PUBLIC :: wood                              ! The Wood function, N = 4.
      50             :   PUBLIC :: chebyquad                         ! The Chebyquad function, 1 <= N.
      51             :   PUBLIC :: leon_cubic_valley                 ! Leon''s cubic valley function, N = 2.
      52             :   PUBLIC :: gregory_karney_tridia_matrix      ! Gregory and Karney''s Tridiagonal Matrix Function, 1 <= N.
      53             :   PUBLIC :: hilbert                           ! The Hilbert function, 1 <= N.
      54             :   PUBLIC :: de_jong_f1                        ! The De Jong Function F1, N = 3.
      55             :   PUBLIC :: de_jong_f2                        ! The De Jong Function F2, N = 2.
      56             :   PUBLIC :: de_jong_f3                        ! The De Jong Function F3 (discontinuous), N = 5.
      57             :   PUBLIC :: de_jong_f4                        ! The De Jong Function F4 (Gaussian noise), N = 30.
      58             :   PUBLIC :: de_jong_f5                        ! The De Jong Function F5, N = 2.
      59             :   PUBLIC :: schaffer_f6                       ! The Schaffer Function F6, N = 2.
      60             :   PUBLIC :: schaffer_f7                       ! The Schaffer Function F7, N = 2.
      61             :   PUBLIC :: goldstein_price_polynomial        ! The Goldstein Price Polynomial, N = 2.
      62             :   PUBLIC :: branin_rcos                       ! The Branin RCOS Function, N = 2.
      63             :   PUBLIC :: shekel_sqrn5                      ! The Shekel SQRN5 Function, N = 4.
      64             :   PUBLIC :: shekel_sqrn7                      ! The Shekel SQRN7 Function, N = 4.
      65             :   PUBLIC :: shekel_sqrn10                     ! The Shekel SQRN10 Function, N = 4.
      66             :   PUBLIC :: six_hump_camel_back_polynomial    ! The Six-Hump Camel-Back Polynomial, N = 2.
      67             :   PUBLIC :: shubert                           ! The Shubert Function, N = 2.
      68             :   PUBLIC :: stuckman                          ! The Stuckman Function, N = 2.
      69             :   PUBLIC :: easom                             ! The Easom Function, N = 2.
      70             :   PUBLIC :: bohachevsky1                      ! The Bohachevsky Function #1, N = 2.
      71             :   PUBLIC :: bohachevsky2                      ! The Bohachevsky Function #2, N = 2.
      72             :   PUBLIC :: bohachevsky3                      ! The Bohachevsky Function #3, N = 2.
      73             :   PUBLIC :: colville_polynomial               ! The Colville Polynomial, N = 4.
      74             :   PUBLIC :: powell3d                          ! The Powell 3D function, N = 3.
      75             :   PUBLIC :: himmelblau                        ! The Himmelblau function, N = 2.
      76             :   ! ------------------------------------------------------------------
      77             :   ! Miscellaneous functions
      78             :   PUBLIC :: griewank                          ! Griewank function, N = 2 or N = 10.
      79             :   PUBLIC :: rosenbrock                        ! Rosenbrock parabolic valley function, N = 2.
      80             :   PUBLIC :: sphere_model                      ! Sphere model, N >= 1.
      81             :   PUBLIC :: rastrigin                         ! Rastrigin function, N >= 2.
      82             :   PUBLIC :: schwefel                          ! Schwefel function, N >= 2.
      83             :   PUBLIC :: ackley                            ! Ackley function, N >= 2.
      84             :   PUBLIC :: michalewicz                       ! Michalewicz function, N >= 2.
      85             :   PUBLIC :: booth                             ! Booth function, N = 2.
      86             :   PUBLIC :: hump                              ! Hump function, N = 2.
      87             :   PUBLIC :: levy                              ! Levy function, N >= 2.
      88             :   PUBLIC :: matyas                            ! Matyas function, N = 2.
      89             :   PUBLIC :: perm                              ! Perm function, N = 4.
      90             :   PUBLIC :: power_sum                         ! Power sum function, N = 4.
      91             :   ! ------------------------------------------------------------------
      92             :   ! test_optimization package of John Burkardt - inputs are x(n) and output f(m), e.g. compare
      93             :   ! rosenbrock = 100.0_dp * (x(2)-x(1)**2)**2 + (1.0_dp-x(1))**2
      94             :   ! rosenbrock_2d(j) = sum((1.0_dp-x(1:m,j))**2) + sum((x(2:m,j)-x(1:m-1,j))**2)
      95             :   PUBLIC :: sphere_model_2d                   !  The sphere model, (M,N).
      96             :   PUBLIC :: axis_parallel_hyper_ellips_2d     !  The axis-parallel hyper-ellipsoid function, (M,N).
      97             :   PUBLIC :: rotated_hyper_ellipsoid_2d        !  The rotated hyper-ellipsoid function, (M,N).
      98             :   PUBLIC :: rosenbrock_2d                     !  Rosenbrock''s valley, (M,N).
      99             :   PUBLIC :: rastrigin_2d                      !  Rastrigin''s function, (M,N).
     100             :   PUBLIC :: schwefel_2d                       !  Schwefel''s function, (M,N).
     101             :   PUBLIC :: griewank_2d                       !  Griewank''s function, (M,N).
     102             :   PUBLIC :: power_sum_2d                      !  The power sum function, (M,N).
     103             :   PUBLIC :: ackley_2d                         !  Ackley''s function, (M,N).
     104             :   PUBLIC :: michalewicz_2d                    !  Michalewicz''s function, (M,N).
     105             :   PUBLIC :: drop_wave_2d                      !  The drop wave function, (M,N).
     106             :   PUBLIC :: deceptive_2d                      !  The deceptive function, (M,N).
     107             :   ! ------------------------------------------------------------------
     108             :   ! test_optimization functions of Kalyanmoy Deb
     109             :   ! found in Deb et al. (2002), Zitzler et al. (2000) and in Matlab Central file exchange
     110             :   ! http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/
     111             :   ! content/TP_NSGA2/
     112             :   public :: dtlz2_3d   !  3-d objective function (spherical                            pareto front)
     113             :   public :: dtlz2_5d   !  5-d objective function (spherical                            pareto front)
     114             :   public :: dtlz2_10d  ! 10-d objective function (spherical                            pareto front)
     115             :   public :: fon_2d     !  2-d objective function (nonconvex                            pareto front)
     116             :   public :: kur_2d     !  2-d objective function (nonconvex,              disconnected pareto front)
     117             :   public :: pol_2d     !  2-d objective function (nonconvex,              disconnected pareto front)
     118             :   public :: sch_2d     !  2-d objective function (   convex                            pareto front)
     119             :   public :: zdt1_2d    !  2-d objective function (   convex                            pareto front)
     120             :   public :: zdt2_2d    !  2-d objective function (nonconvex                            pareto front)
     121             :   public :: zdt3_2d    !  2-d objective function (   convex,              disconnected pareto front)
     122             :   public :: zdt4_2d    !  2-d objective function (nonconvex                            pareto front)
     123             :   public :: zdt6_2d    !  2-d objective function (nonconvex, nonuniformly disconnected pareto front)
     124             : 
     125             :   ! routines added to be compatible with testing framework
     126             :   public :: ackley_objective, griewank_objective, eval_dummy
     127             : 
     128             : CONTAINS
     129             : 
     130             :   ! ------------------------------------------------------------------
     131             :   !
     132             :   !  Simple quadratic, (x-2)^2+1
     133             :   !  Solution: x = 2.0
     134             :   !  With Brent method:
     135             :   !   A,  X*,  B:   1.9999996       2.0000000       2.0000004
     136             :   !  FA, FX*, FB:   1.0000000       1.0000000       1.0000000
     137             :   !
     138             :   !  Licensing:
     139             :   !
     140             :   !    This code is distributed under the GNU LGPL license.
     141             :   !
     142             :   !  Modified:
     143             :   !
     144             :   !    25 February 2002
     145             :   !
     146             :   !  Author:
     147             :   !
     148             :   !    John Burkardt
     149             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     150             :   !
     151             :   !  Parameters:
     152             :   !
     153             :   !    Input, real(dp) X, the argument of the objective function.
     154             :   !
     155             : 
     156           0 :   function quadratic( x )
     157             : 
     158             :     implicit none
     159             : 
     160             :     real(dp), dimension(:), intent(in) :: x
     161             :     real(dp)                           :: quadratic
     162             : 
     163           0 :     if (size(x,1) .gt. 1_i4) stop 'quadratic: Input has to be array of size 1'
     164           0 :     quadratic = ( x(1) - 2.0_dp ) * ( x(1) - 2.0_dp ) + 1.0_dp
     165             : 
     166           0 :   end function quadratic
     167             : 
     168             :   ! ------------------------------------------------------------------
     169             :   !
     170             :   !  Quadratic plus exponential, x^2 + e^(-x)
     171             :   !  Solution: x = 0.35173370
     172             :   !  With Brent method:
     173             :   !   A,  X*,  B:  0.35173337      0.35173370      0.35173404
     174             :   !  FA, FX*, FB:  0.82718403      0.82718403      0.82718403
     175             :   !
     176             :   !  Licensing:
     177             :   !
     178             :   !    This code is distributed under the GNU LGPL license.
     179             :   !
     180             :   !  Modified:
     181             :   !
     182             :   !    26 February 2002
     183             :   !
     184             :   !  Author:
     185             :   !
     186             :   !    John Burkardt
     187             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     188             :   !
     189             :   !  Reference:
     190             :   !
     191             :   !    LE Scales,
     192             :   !    Introduction to Non-Linear Optimization,
     193             :   !    Springer, 1985.
     194             :   !
     195             :   !  Parameters:
     196             :   !
     197             :   !    Input, real(dp) X, the argument of the objective function.
     198             :   !
     199             : 
     200           0 :   function quadratic_exponential( x )
     201             : 
     202             :     implicit none
     203             : 
     204             :     real(dp), dimension(:), intent(in) :: x
     205             :     real(dp)                           :: quadratic_exponential
     206             : 
     207           0 :     if (size(x,1) .gt. 1_i4) stop 'quadratic_exponential: Input has to be array of size 1'
     208           0 :     quadratic_exponential = x(1) * x(1) + exp ( - x(1) )
     209             : 
     210           0 :   end function quadratic_exponential
     211             : 
     212             :   ! ------------------------------------------------------------------
     213             :   !
     214             :   !  Quartic, x^4 + 2x^2 + x + 3
     215             :   !  Solution: x = -0.23673291
     216             :   !  With Brent method:
     217             :   !   A,  X*,  B: -0.23673324     -0.23673291     -0.23673257
     218             :   !  FA, FX*, FB:   2.8784928       2.8784928       2.8784928
     219             :   !
     220             :   !  Licensing:
     221             :   !
     222             :   !    This code is distributed under the GNU LGPL license.
     223             :   !
     224             :   !  Modified:
     225             :   !
     226             :   !    26 February 2002
     227             :   !
     228             :   !  Author:
     229             :   !
     230             :   !    John Burkardt
     231             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     232             :   !
     233             :   !  Reference:
     234             :   !
     235             :   !    LE Scales,
     236             :   !    Introduction to Non-Linear Optimization,
     237             :   !    Springer, 1985.
     238             :   !
     239             :   !  Parameters:
     240             :   !
     241             :   !    Input, real(dp) X, the argument of the objective function.
     242             :   !
     243             : 
     244           0 :   function quartic( x )
     245             : 
     246             :     implicit none
     247             : 
     248             :     real(dp), dimension(:), intent(in) :: x
     249             :     real(dp)                           :: quartic
     250             : 
     251           0 :     if (size(x,1) .gt. 1_i4) stop 'quartic: Input has to be array of size 1'
     252           0 :     quartic = ( ( x(1) * x(1) + 2.0_dp ) * x(1) + 1.0_dp ) * x(1) + 3.0_dp
     253             : 
     254           0 :   end function quartic
     255             : 
     256             :   ! ------------------------------------------------------------------
     257             :   !
     258             :   !  Steep valley, e^x + 1/(100x)
     259             :   !  Solution:
     260             :   !       if x >  0.0 :   x = 0.95344636E-01
     261             :   !       if x < -0.1 :   x = -8.99951
     262             :   !  Search domain: x <= -0.1
     263             :   !
     264             :   !  With Brent method:
     265             :   !   A,  X*,  B:  0.95344301E-01  0.95344636E-01  0.95344971E-01
     266             :   !  FA, FX*, FB:   1.2049206       1.2049206       1.2049206
     267             :   !
     268             :   !  Discussion:
     269             :   !
     270             :   !    This function has a pole at x = 0,
     271             :   !    near which
     272             :   !       f -> negative infinity for x -> 0-0 (left) and
     273             :   !       f -> positive infinity for x -> 0+0 (right)
     274             :   !    f has a local maximum at x ~ -0.105412 .
     275             :   !
     276             :   !  Licensing:
     277             :   !
     278             :   !    This code is distributed under the GNU LGPL license.
     279             :   !
     280             :   !  Modified:
     281             :   !
     282             :   !    26 February 2002
     283             :   !
     284             :   !  Author:
     285             :   !
     286             :   !    John Burkardt
     287             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     288             :   !
     289             :   !  Reference:
     290             :   !
     291             :   !    LE Scales,
     292             :   !    Introduction to Non-Linear Optimization,
     293             :   !    Springer, 1985.
     294             :   !
     295             :   !  Parameters:
     296             :   !
     297             :   !    Input, real(dp) X, the argument of the objective function.
     298             :   !
     299             : 
     300           0 :   function steep_valley( x )
     301             : 
     302             :     implicit none
     303             : 
     304             :     real(dp), dimension(:), intent(in) :: x
     305             :     real(dp)                           :: steep_valley
     306             : 
     307           0 :     if (size(x,1) .gt. 1_i4) stop 'steep_valley: Input has to be array of size 1'
     308           0 :     steep_valley = exp ( x(1) ) + 0.01_dp / x(1)
     309             : 
     310           0 :     steep_valley = steep_valley + 0.0009877013_dp
     311             : 
     312           0 :   end function steep_valley
     313             : 
     314             :   ! ------------------------------------------------------------------
     315             :   !
     316             :   !  Steep valley2, e^x - 2x + 1/(100x) - 1/(1000000x^2)
     317             :   !
     318             :   !  Solution: x = 0.70320487
     319             :   !  Search domain: 0.001 <= x
     320             :   !
     321             :   !  With Brent method:
     322             :   !   A,  X*,  B:  0.70320453      0.70320487      0.70320521
     323             :   !  FA, FX*, FB:  0.62802572      0.62802572      0.62802572
     324             :   !
     325             :   !  Licensing:
     326             :   !
     327             :   !    This code is distributed under the GNU LGPL license.
     328             :   !
     329             :   !  Modified:
     330             :   !
     331             :   !    26 February 2002
     332             :   !
     333             :   !  Author:
     334             :   !
     335             :   !    John Burkardt
     336             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     337             :   !
     338             :   !  Reference:
     339             :   !
     340             :   !    LE Scales,
     341             :   !    Introduction to Non-Linear Optimization,
     342             :   !    Springer, 1985.
     343             :   !
     344             :   !  Parameters:
     345             :   !
     346             :   !    Input, real(dp) X, the argument of the objective function.
     347             :   !
     348             : 
     349           0 :   function steep_valley2( x )
     350             : 
     351             :     implicit none
     352             : 
     353             :     real(dp), dimension(:), intent(in) :: x
     354             :     real(dp)                           :: steep_valley2
     355             : 
     356           0 :     if (size(x,1) .gt. 1_i4) stop 'steep_valley2: Input has to be array of size 1'
     357           0 :     steep_valley2 = exp( x(1) ) - 2.0_dp * x(1) + 0.01_dp / x(1) - 0.000001_dp / x(1) / x(1)
     358             : 
     359           0 :   end function steep_valley2
     360             : 
     361             :   ! ------------------------------------------------------------------
     362             :   !
     363             :   !  The dying snake, ( x + sin(x) ) * e^(-x^2)
     364             :   !  Solution: x = -0.67957876
     365             :   !  With Brent method:
     366             :   !   A,  X*,  B: -0.67957911     -0.67957876     -0.67957842
     367             :   !  FA, FX*, FB: -0.82423940     -0.82423940     -0.82423940
     368             :   !
     369             :   !  Licensing:
     370             :   !
     371             :   !    This code is distributed under the GNU LGPL license.
     372             :   !
     373             :   !  Modified:
     374             :   !
     375             :   !    26 February 2002
     376             :   !
     377             :   !  Author:
     378             :   !
     379             :   !    John Burkardt
     380             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     381             :   !
     382             :   !  Reference:
     383             :   !
     384             :   !    Richard Brent,
     385             :   !    Algorithms for Minimization Without Derivatives,
     386             :   !    Prentice Hall 1973,
     387             :   !    Reprinted Dover, 2002
     388             :   !
     389             :   !  Parameters:
     390             :   !
     391             :   !    Input, real(dp) X, the argument of the objective function.
     392             :   !
     393             : 
     394           0 :   function dying_snake(x)
     395             : 
     396             :     implicit none
     397             : 
     398             :     real(dp), dimension(:), intent(in) :: x
     399             :     real(dp)                           :: dying_snake
     400             : 
     401           0 :     if (size(x,1) .gt. 1_i4) stop 'dying_snake: Input has to be array of size 1'
     402           0 :     dying_snake = ( x(1) + sin ( x(1) ) ) * exp ( - x(1) * x(1) )
     403             : 
     404           0 :     dying_snake = dying_snake + 0.8242393985_dp
     405             : 
     406           0 :   end function dying_snake
     407             : 
     408             :   ! ------------------------------------------------------------------
     409             :   !
     410             :   !  The "Thin Pole", 3x^2+1+log((pi-x)^2)/pi^4
     411             :   !  Solution:
     412             :   !     x    = 0.00108963
     413             :   !     f(x) = 1.0235
     414             :   !
     415             :   !  With Brent method:
     416             :   !   A,  X*,  B:   2.0000000       2.0000007       2.0000011
     417             :   !  FA, FX*, FB:   13.002719       13.002727       13.002732
     418             :   !
     419             :   !  Discussion:
     420             :   !
     421             :   !    This function looks positive, but has a pole at x = pi,
     422             :   !    near which f -> negative infinity, and has two zeroes nearby.
     423             :   !    None of this will show up computationally.
     424             :   !
     425             :   !  Licensing:
     426             :   !
     427             :   !    This code is distributed under the GNU LGPL license.
     428             :   !
     429             :   !  Modified:
     430             :   !
     431             :   !    19 February 2003
     432             :   !
     433             :   !  Author:
     434             :   !
     435             :   !    John Burkardt
     436             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     437             :   !
     438             :   !  Reference:
     439             :   !
     440             :   !    Arnold Krommer, Christoph Ueberhuber,
     441             :   !    Numerical Integration on Advanced Systems,
     442             :   !    Springer, 1994, pages 185-186.
     443             :   !
     444             :   !  Parameters:
     445             :   !
     446             :   !    Input, real(dp) X, the argument of the objective function.
     447             :   !
     448             : 
     449           0 :   function thin_pole(x)
     450             : 
     451           0 :     use mo_constants, only: pi_dp
     452             :     use mo_utils,     only: eq
     453             : 
     454             :     implicit none
     455             : 
     456             :     real(dp), dimension(:), intent(in) :: x
     457             :     real(dp)                           :: thin_pole
     458             : 
     459           0 :     if (size(x,1) .gt. 1_i4) stop 'thin_pole: Input has to be array of size 1'
     460           0 :     if ( eq(x(1),pi_dp) ) then
     461             :        thin_pole = - 10000.0_dp
     462             :     else
     463           0 :        thin_pole = 3.0_dp * x(1) * x(1) + 1.0_dp + ( log ( ( x(1) - pi_dp ) * ( x(1) - pi_dp ) ) ) / pi_dp**4
     464             :     end if
     465             : 
     466           0 :   end function thin_pole
     467             : 
     468             :   ! ------------------------------------------------------------------
     469             :   !
     470             :   !  The oscillatory parabola x^2 - 10*sin(x^2-3x+2)
     471             :   !  Solution:
     472             :   !     x    = 0.146623
     473             :   !     f(x) = -9.97791
     474             :   !  With Brent method:
     475             :   !   A,  X*,  B:  -1.3384524      -1.3384521      -1.3384517
     476             :   !  FA, FX*, FB:  -8.1974224      -8.1974224      -8.1974224
     477             :   !
     478             :   !  Discussion:
     479             :   !
     480             :   !    This function is oscillatory, with many local minima.
     481             :   !
     482             :   !  Licensing:
     483             :   !
     484             :   !    This code is distributed under the GNU LGPL license.
     485             :   !
     486             :   !  Modified:
     487             :   !
     488             :   !    25 January 2008
     489             :   !
     490             :   !  Author:
     491             :   !
     492             :   !    John Burkardt
     493             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     494             :   !
     495             :   !  Parameters:
     496             :   !
     497             :   !    Input, real(dp) X, the argument of the objective function.
     498             :   !
     499             : 
     500           0 :   function oscillatory_parabola(x)
     501             : 
     502             :     implicit none
     503             : 
     504             :     real(dp), dimension(:), intent(in) :: x
     505             :     real(dp)                           :: oscillatory_parabola
     506             : 
     507           0 :     if (size(x,1) .gt. 1_i4) stop 'oscillatory_parabola: Input has to be array of size 1'
     508           0 :     oscillatory_parabola = x(1) * x(1) - 10.0_dp * sin ( x(1) * x(1) - 3.0_dp * x(1) + 2.0_dp )
     509             : 
     510           0 :     oscillatory_parabola = oscillatory_parabola + 9.9779149346_dp
     511             : 
     512           0 :   end function oscillatory_parabola
     513             : 
     514             :   ! ------------------------------------------------------------------
     515             :   !
     516             :   !  The cosine combo cos(x)+5cos(1.6x)-2cos(2x)+5cos(4.5x)+7cos(9x)
     517             :   !  Solution:
     518             :   !     x    = -21.9443 + 62.831853 * k   , k = Integer
     519             :   !     x    =  21.9443 - 62.831853 * k   , k = Integer
     520             :   !     f(x) = -14.6772
     521             :   !
     522             :   !  With Brent method:
     523             :   !   A,  X*,  B:   1.0167817       1.0167821       1.0167824
     524             :   !  FA, FX*, FB:  -6.2827509      -6.2827509      -6.2827509
     525             :   !
     526             :   !  Discussion:
     527             :   !
     528             :   !    This function is symmetric, oscillatory, and has many local minima.
     529             :   !
     530             :   !    The function has a local minimum at 1.0167817.
     531             :   !
     532             :   !    The global optimum which function value -14.6772
     533             :   !    appears infinite times.
     534             :   !
     535             :   !  Licensing:
     536             :   !
     537             :   !    This code is distributed under the GNU LGPL license.
     538             :   !
     539             :   !  Modified:
     540             :   !
     541             :   !    09 February 2009
     542             :   !
     543             :   !  Author:
     544             :   !
     545             :   !    John Burkardt
     546             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     547             :   !
     548             :   !  Reference:
     549             :   !
     550             :   !    Isabel Beichl, Dianne O'Leary, Francis Sullivan,
     551             :   !    Monte Carlo Minimization and Counting: One, Two, Too Many,
     552             :   !    Computing in Science and Engineering,
     553             :   !    Volume 9, Number 1, January/February 2007.
     554             :   !
     555             :   !    Dianne O'Leary,
     556             :   !    Scientific Computing with Case Studies,
     557             :   !    SIAM, 2008,
     558             :   !    ISBN13: 978-0-898716-66-5,
     559             :   !    LC: QA401.O44.
     560             :   !
     561             :   !  Parameters:
     562             :   !
     563             :   !    Input, real(dp) X, the argument of the objective function.
     564             :   !
     565             : 
     566           0 :   function cosine_combo(x)
     567             : 
     568             :     implicit none
     569             : 
     570             :     real(dp), dimension(:), intent(in) :: x
     571             :     real(dp)                           :: cosine_combo
     572             : 
     573           0 :     if (size(x,1) .gt. 1_i4) stop 'cosine_combo: Input has to be array of size 1'
     574           0 :     cosine_combo =   cos (         x(1) ) &
     575             :          + 5.0_dp * cos ( 1.6_dp * x(1) ) &
     576             :          - 2.0_dp * cos ( 2.0_dp * x(1) ) &
     577             :          + 5.0_dp * cos ( 4.5_dp * x(1) ) &
     578           0 :          + 7.0_dp * cos ( 9.0_dp * x(1) )
     579             : 
     580           0 :     cosine_combo = cosine_combo + 14.6771885214_dp
     581             : 
     582           0 :   end function cosine_combo
     583             : 
     584             :   ! ------------------------------------------------------------------
     585             :   !
     586             :   !  abs1, 1 + |3x-1|
     587             :   !  Solution: x = 1./3.
     588             :   !  With Brent method:
     589             :   !   A,  X*,  B:  0.33333299      0.33333351      0.33333385
     590             :   !  FA, FX*, FB:   1.0000010       1.0000005       1.0000015
     591             :   !
     592             :   !  Licensing:
     593             :   !
     594             :   !    This code is distributed under the GNU LGPL license.
     595             :   !
     596             :   !  Modified:
     597             :   !
     598             :   !    03 February 2012
     599             :   !
     600             :   !  Author:
     601             :   !
     602             :   !    John Burkardt
     603             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     604             :   !
     605             :   !  Parameters:
     606             :   !
     607             :   !    Input, real(dp) X, the argument of the objective function.
     608             :   !
     609             : 
     610           0 :   function abs1(x)
     611             : 
     612             :     implicit none
     613             : 
     614             :     real(dp), dimension(:), intent(in) :: x
     615             :     real(dp)                           :: abs1
     616             : 
     617           0 :     if (size(x,1) .gt. 1_i4) stop 'abs1: Input has to be array of size 1'
     618           0 :     abs1 = 1.0_dp + abs ( 3.0_dp * x(1) - 1.0_dp )
     619             : 
     620           0 :   end function abs1
     621             : 
     622             :   ! ------------------------------------------------------------------
     623             :   !
     624             :   !  R8_AINT truncates an R8 argument to an integer.
     625             :   !
     626             :   !  Licensing:
     627             :   !
     628             :   !    This code is distributed under the GNU LGPL license.
     629             :   !
     630             :   !  Modified:
     631             :   !
     632             :   !    18 October 2011
     633             :   !
     634             :   !  Author:
     635             :   !
     636             :   !    John Burkardt.
     637             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     638             :   !
     639             :   !  Parameters:
     640             :   !
     641             :   !    Input, real(dp) X, the argument.
     642             :   !
     643             :   !    Output, real(dp) VALUE, the truncated version of X.
     644             :   !
     645             : 
     646           0 :   function r8_aint( x )
     647             : 
     648             :     implicit none
     649             : 
     650             :     real(dp) :: r8_aint
     651           0 :     real(dp) :: val
     652             :     real(dp) :: x
     653             : 
     654           0 :     if ( x < 0.0_dp ) then
     655           0 :        val = -int( abs ( x ) )
     656             :     else
     657           0 :        val =  int( abs ( x ) )
     658             :     end if
     659             : 
     660           0 :     r8_aint = val
     661             : 
     662           0 :   end function r8_aint
     663             : 
     664             :   !*****************************************************************************80
     665             :   !
     666             :   !! NORMAL_01_SAMPLE samples the standard Normal PDF.
     667             :   !
     668             :   !  Discussion:
     669             :   !
     670             :   !    The standard normal distribution has mean 0 and standard
     671             :   !    deviation 1.
     672             :   !
     673             :   !  Method:
     674             :   !
     675             :   !    The Box-Muller method is used.
     676             :   !
     677             :   !  Licensing:
     678             :   !
     679             :   !    This code is distributed under the GNU LGPL license.
     680             :   !
     681             :   !  Modified:
     682             :   !
     683             :   !    01 December 2000
     684             :   !
     685             :   !  Author:
     686             :   !
     687             :   !    John Burkardt
     688             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     689             :   !
     690             :   !  Parameters:
     691             :   !
     692             :   !    Output, real(dp) X, a sample of the PDF.
     693             :   !
     694             : 
     695           0 :   subroutine normal_01_sample ( x )
     696             : 
     697           0 :     use mo_constants, only: pi_dp
     698             :     use mo_utils, only: le
     699             : 
     700             :     implicit none
     701             : 
     702             :     integer(i4), save :: iset = -1
     703           0 :     real(dp) v1
     704           0 :     real(dp) v2
     705             :     real(dp) x
     706             :     real(dp), save :: xsave = 0.0_dp
     707             : 
     708           0 :     if ( iset == -1 ) then
     709           0 :        call random_seed ( )
     710           0 :        iset = 0
     711             :     end if
     712             : 
     713           0 :     if ( iset == 0 ) then
     714             : 
     715           0 :        call random_number ( harvest = v1 )
     716             : 
     717           0 :        if ( le(v1,0.0_dp) ) then
     718           0 :           write ( *, '(a)' ) ' '
     719           0 :           write ( *, '(a)' ) 'NORMAL_01_SAMPLE - Fatal error!'
     720           0 :           write ( *, '(a)' ) '  V1 <= 0.'
     721           0 :           write ( *, '(a,g14.6)' ) '  V1 = ', v1
     722           0 :           stop
     723             :        end if
     724             : 
     725           0 :        call random_number ( harvest = v2 )
     726             : 
     727           0 :        if ( le(v2,0.0_dp) ) then
     728           0 :           write ( *, '(a)' ) ' '
     729           0 :           write ( *, '(a)' ) 'NORMAL_01_SAMPLE - Fatal error!'
     730           0 :           write ( *, '(a)' ) '  V2 <= 0.'
     731           0 :           write ( *, '(a,g14.6)' ) '  V2 = ', v2
     732           0 :           stop
     733             :        end if
     734             : 
     735           0 :        x = sqrt ( - 2.0_dp * log ( v1 ) ) * cos ( 2.0_dp * pi_dp * v2 )
     736             : 
     737           0 :        xsave = sqrt ( - 2.0_dp * log ( v1 ) ) * sin ( 2.0_dp * PI_dp * v2 )
     738             : 
     739           0 :        iset = 1
     740             : 
     741             :     else
     742             : 
     743           0 :        x = xsave
     744           0 :        iset = 0
     745             : 
     746             :     end if
     747             : 
     748           0 :     return
     749           0 :   end subroutine normal_01_sample
     750             : 
     751             :   ! ------------------------------------------------------------------
     752             :   !
     753             :   ! The Fletcher-Powell helical valley function, N = 3.
     754             :   ! Solution: x(1:3) = (/ 1.0_dp, 0.0_dp, 0.0_dp /)
     755             :   !
     756             :   !  Licensing:
     757             :   !
     758             :   !    This code is distributed under the GNU LGPL license.
     759             :   !
     760             :   !  Modified:
     761             :   !
     762             :   !    15 March 2000
     763             :   !
     764             :   !  Author:
     765             :   !
     766             :   !    John Burkardt
     767             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     768             :   !
     769             :   !  Reference:
     770             :   !
     771             :   !    Richard Brent,
     772             :   !    Algorithms for Minimization with Derivatives,
     773             :   !    Dover, 2002,
     774             :   !    ISBN: 0-486-41998-3,
     775             :   !    LC: QA402.5.B74.
     776             :   !
     777             :   !  Parameters:
     778             :   !
     779             :   !    Input, real(dp) :: X(N), the argument of the objective function.
     780             :   !
     781             : 
     782           0 :   function fletcher_powell_helical_valley(x)
     783             : 
     784           0 :     use mo_constants, only: pi_dp
     785             : 
     786             :     implicit none
     787             : 
     788             :     !    integer(i4) :: n
     789             : 
     790             :     real(dp) :: fletcher_powell_helical_valley
     791           0 :     real(dp) :: th
     792             :     real(dp), dimension(:), intent(in) :: x
     793             : 
     794           0 :     if ( 0.0_dp < x(1) ) then
     795           0 :        th = 0.5_dp * atan ( x(2) / x(1) ) / pi_dp
     796           0 :     else if ( x(1) < 0.0_dp ) then
     797           0 :        th = 0.5_dp * atan ( x(2) / x(1) ) / pi_dp + 0.5_dp
     798           0 :     else if ( 0.0_dp < x(2) ) then
     799             :        th = 0.25_dp
     800           0 :     else if ( x(2) < 0.0_dp ) then
     801             :        th = - 0.25_dp
     802             :     else
     803           0 :        th = 0.0_dp
     804             :     end if
     805             :     !call p01_th ( x, th )
     806             : 
     807           0 :     fletcher_powell_helical_valley = 100.0_dp * ( x(3) - 10.0_dp * th )**2 &
     808           0 :          + 100.0_dp * ( sqrt ( x(1) * x(1) + x(2) * x(2) ) - 1.0_dp )**2 &
     809           0 :          + x(3) * x(3)
     810             : 
     811           0 :   end function fletcher_powell_helical_valley
     812             : 
     813             :   ! ------------------------------------------------------------------
     814             :   !
     815             :   ! The Biggs EXP6 function, N = 6.
     816             :   ! Solution: x(1:6) = (/ 1.0_dp, 10.0_dp, 1.0_dp, 5.0_dp, 4.0_dp, 3.0_dp /)
     817             :   !           at f(x*) = 0.0
     818             :   !
     819             :   !  Licensing:
     820             :   !
     821             :   !    This code is distributed under the GNU LGPL license.
     822             :   !
     823             :   !  Modified:
     824             :   !
     825             :   !    04 May 2000
     826             :   !
     827             :   !  Author:
     828             :   !
     829             :   !    John Burkardt
     830             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     831             :   !
     832             :   !  Parameters:
     833             :   !
     834             :   !    Input, real(dp) :: X(N), the argument of the objective function.
     835             :   !
     836             : 
     837           0 :   function biggs_exp6(x)
     838             : 
     839             :     implicit none
     840             : 
     841             :     !    integer(i4) :: n
     842             : 
     843           0 :     real(dp) :: c
     844             :     real(dp) :: biggs_exp6
     845           0 :     real(dp) :: fi
     846             :     integer(i4) ::i
     847             :     real(dp), dimension(:), intent(in) :: x
     848             : 
     849           0 :     biggs_exp6 = 0.0_dp
     850             : 
     851           0 :     do i = 1, 13
     852             : 
     853           0 :        c = - real ( i, dp ) / 10.0_dp
     854             : 
     855           0 :        fi = x(3)     * exp ( c * x(1) )         - x(4) * exp ( c * x(2) ) &
     856           0 :             + x(6)     * exp ( c * x(5) )         -        exp ( c ) &
     857           0 :             + 5.0_dp  * exp ( 10.0_dp * c ) - 3.0_dp  * exp ( 4.0_dp * c )
     858             : 
     859           0 :        biggs_exp6 = biggs_exp6 + fi * fi
     860             : 
     861             :     end do
     862             : 
     863           0 :   end function biggs_exp6
     864             : 
     865             :   ! ------------------------------------------------------------------
     866             :   !
     867             :   ! The Gaussian function, N = 3.
     868             :   ! Search domain: -Pi <= xi <= Pi, i = 1, 2, 3.
     869             :   ! Solution:
     870             :   !     x(1:n) = (/ 0.39895613783875655_dp, 1.0000190844878036_dp, 0.0_dp /)
     871             :   !     at f(x*) = 0.0
     872             :   !     found with Mathematica
     873             :   !
     874             :   !  Licensing:
     875             :   !
     876             :   !    This code is distributed under the GNU LGPL license.
     877             :   !
     878             :   !  Modified:
     879             :   !
     880             :   !    24 March 2000
     881             :   !
     882             :   !  Author:
     883             :   !
     884             :   !    John Burkardt
     885             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     886             :   !
     887             :   !  Parameters:
     888             :   !
     889             :   !    Input, real(dp) :: X(N), the argument of the objective function.
     890             :   !
     891             : 
     892           0 :   function gaussian(x)
     893             : 
     894             :     implicit none
     895             : 
     896             :     !    integer(i4) :: n
     897             : 
     898             :     real(dp), dimension(:), intent(in) :: x
     899             :     real(dp)                           :: gaussian
     900             : 
     901             :     integer(i4) :: i
     902           0 :     real(dp)    :: t
     903           0 :     real(dp)    :: y(15)
     904             : 
     905             :     y(1:15) = (/ 0.0009_dp, 0.0044_dp, 0.0175_dp, 0.0540_dp, 0.1295_dp, &
     906             :          0.2420_dp, 0.3521_dp, 0.3989_dp, 0.3521_dp, 0.2420_dp, &
     907           0 :          0.1295_dp, 0.0540_dp, 0.0175_dp, 0.0044_dp, 0.0009_dp /)
     908             : 
     909           0 :     gaussian = 0.0_dp
     910             : 
     911           0 :     do i = 1, 15
     912             : 
     913             :        ! avoiding underflow
     914           0 :        t =  - 0.5_dp * x(2) * &
     915           0 :             ( 3.5_dp - 0.5_dp * real ( i - 1, dp ) - x(3) )**2
     916           0 :        if ( t .lt. -709._dp ) then
     917           0 :           t = -y(i)
     918             :        else
     919           0 :           t = x(1) * exp ( t ) - y(i)
     920             :        end if
     921             : 
     922           0 :        gaussian = gaussian + t * t
     923             : 
     924             :     end do
     925             : 
     926           0 :   end function gaussian
     927             : 
     928             :   ! ------------------------------------------------------------------
     929             :   !
     930             :   ! The Powell badly scaled function, N = 2.
     931             :   ! Solution: x(1:2) = (/ 1.098159D-05, 9.106146_dp /)
     932             :   !
     933             :   !  Licensing:
     934             :   !
     935             :   !    This code is distributed under the GNU LGPL license.
     936             :   !
     937             :   !  Modified:
     938             :   !
     939             :   !    04 May 2000
     940             :   !
     941             :   !  Author:
     942             :   !
     943             :   !    John Burkardt
     944             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     945             :   !
     946             :   !  Reference:
     947             :   !
     948             :   !    Richard Brent,
     949             :   !    Algorithms for Minimization with Derivatives,
     950             :   !    Dover, 2002,
     951             :   !    ISBN: 0-486-41998-3,
     952             :   !    LC: QA402.5.B74.
     953             :   !
     954             :   !  Parameters:
     955             :   !
     956             :   !    Input, real(dp) :: X(N), the argument of the objective function.
     957             :   !
     958             : 
     959           0 :   function powell_badly_scaled(x)
     960             : 
     961             :     implicit none
     962             : 
     963             :     !    integer(i4) :: n
     964             : 
     965             :     real(dp) :: powell_badly_scaled
     966           0 :     real(dp) :: f1
     967           0 :     real(dp) :: f2
     968             :     real(dp), dimension(:), intent(in) :: x
     969             : 
     970           0 :     f1 = 10000.0_dp * x(1) * x(2) - 1.0_dp
     971           0 :     f2 = exp ( - x(1) ) + exp ( - x(2) ) - 1.0001_dp
     972             : 
     973           0 :     powell_badly_scaled = f1 * f1 + f2 * f2
     974             : 
     975           0 :   end function powell_badly_scaled
     976             : 
     977             :   ! ------------------------------------------------------------------
     978             :   !
     979             :   ! The Box 3-dimensional function, N = 3.
     980             :   ! Solution: x(1:3) = (/ 1.0_dp, 10.0_dp, 1.0_dp /)
     981             :   ! seems to be not the only solution
     982             :   !
     983             :   !  Discussion:
     984             :   !
     985             :   !    The function is formed by the sum of squares of 10 separate terms.
     986             :   !
     987             :   !  Licensing:
     988             :   !
     989             :   !    This code is distributed under the GNU LGPL license.
     990             :   !
     991             :   !  Modified:
     992             :   !
     993             :   !    04 May 2000
     994             :   !
     995             :   !  Author:
     996             :   !
     997             :   !    John Burkardt
     998             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
     999             :   !
    1000             :   !  Parameters:
    1001             :   !
    1002             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1003             :   !
    1004             : 
    1005           0 :   function box_3dimensional(x)
    1006             : 
    1007             :     implicit none
    1008             : 
    1009             :     !    integer(i4) :: n
    1010             : 
    1011           0 :     real(dp) :: c
    1012             :     real(dp) :: box_3dimensional
    1013           0 :     real(dp) :: fi
    1014             :     integer(i4) ::i
    1015             :     real(dp), dimension(:), intent(in) :: x
    1016             : 
    1017           0 :     box_3dimensional = 0.0_dp
    1018             : 
    1019           0 :     do i = 1, 10
    1020             : 
    1021           0 :        c = - real ( i, dp ) / 10.0_dp
    1022             : 
    1023           0 :        fi = exp ( c * x(1) ) - exp ( c * x(2) ) - x(3) * &
    1024           0 :             ( exp ( c ) - exp ( 10.0_dp * c ) )
    1025             : 
    1026           0 :        box_3dimensional = box_3dimensional + fi * fi
    1027             : 
    1028             :     end do
    1029             : 
    1030           0 :   end function box_3dimensional
    1031             : 
    1032             :   ! ------------------------------------------------------------------
    1033             :   !
    1034             :   ! The variably dimensioned function, 1 <= N.
    1035             :   ! Solution: x(1:n) = 1.0_dp
    1036             :   !
    1037             :   !  Licensing:
    1038             :   !
    1039             :   !    This code is distributed under the GNU LGPL license.
    1040             :   !
    1041             :   !  Modified:
    1042             :   !
    1043             :   !    05 May 2000
    1044             :   !
    1045             :   !  Author:
    1046             :   !
    1047             :   !    John Burkardt
    1048             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1049             :   !
    1050             :   !  Reference:
    1051             :   !
    1052             :   !    Richard Brent,
    1053             :   !    Algorithms for Minimization with Derivatives,
    1054             :   !    Dover, 2002,
    1055             :   !    ISBN: 0-486-41998-3,
    1056             :   !    LC: QA402.5.B74.
    1057             :   !
    1058             :   !  Parameters:
    1059             :   !
    1060             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1061             :   !
    1062             : 
    1063           0 :   function variably_dimensioned(x)
    1064             : 
    1065             :     implicit none
    1066             : 
    1067             :     integer(i4) :: n
    1068             : 
    1069             :     real(dp) :: variably_dimensioned
    1070           0 :     real(dp) :: f1
    1071           0 :     real(dp) :: f2
    1072             :     integer(i4) ::i
    1073             :     real(dp), dimension(:), intent(in) :: x
    1074             : 
    1075           0 :     n = size(x)
    1076           0 :     f1 = 0.0_dp
    1077           0 :     do i = 1, n
    1078           0 :        f1 = f1 + real ( i, dp ) * ( x(i) - 1.0_dp )
    1079             :     end do
    1080             : 
    1081             :     f2 = 0.0_dp
    1082           0 :     do i = 1, n
    1083           0 :        f2 = f2 + ( x(i) - 1.0_dp )**2
    1084             :     end do
    1085             : 
    1086           0 :     variably_dimensioned = f1 * f1 * ( 1.0_dp + f1 * f1 ) + f2
    1087             : 
    1088           0 :   end function variably_dimensioned
    1089             : 
    1090             :   ! ------------------------------------------------------------------
    1091             :   !
    1092             :   ! The Watson function, 2 <= N.
    1093             :   ! Solution: n==6: x(1:n) = (/ -0.015725_dp, 1.012435_dp, -0.232992_dp,
    1094             :   !                             1.260430_dp, -1.513729_dp, 0.992996_dp /)
    1095             :   !           n==9  x(1:n) = (/ -0.000015_dp, 0.999790_dp, 0.014764_dp,
    1096             :   !                              0.146342_dp, 1.000821_dp, -2.617731_dp,
    1097             :   !                              4.104403_dp, -3.143612_dp, 1.052627_dp /)
    1098             :   !           else unknown
    1099             :   !
    1100             :   !  Discussion:
    1101             :   !
    1102             :   !    For N = 9, the problem of minimizing the Watson function is
    1103             :   !    very ill conditioned.
    1104             :   !
    1105             :   !  Licensing:
    1106             :   !
    1107             :   !    This code is distributed under the GNU LGPL license.
    1108             :   !
    1109             :   !  Modified:
    1110             :   !
    1111             :   !    15 March 2000
    1112             :   !
    1113             :   !  Author:
    1114             :   !
    1115             :   !    John Burkardt
    1116             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1117             :   !
    1118             :   !  Reference:
    1119             :   !
    1120             :   !    Richard Brent,
    1121             :   !    Algorithms for Minimization with Derivatives,
    1122             :   !    Dover, 2002,
    1123             :   !    ISBN: 0-486-41998-3,
    1124             :   !    LC: QA402.5.B74.
    1125             :   !
    1126             :   !  Parameters:
    1127             :   !
    1128             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1129             :   !
    1130             : 
    1131           0 :   function watson(x)
    1132             : 
    1133             :     implicit none
    1134             : 
    1135             :     integer(i4) :: n
    1136             : 
    1137           0 :     real(dp) :: d
    1138             :     real(dp) :: watson
    1139             :     integer(i4) ::i
    1140             :     integer(i4) ::j
    1141           0 :     real(dp) :: s1
    1142           0 :     real(dp) :: s2
    1143             :     real(dp), dimension(:), intent(in) :: x
    1144             : 
    1145           0 :     n = size(x)
    1146           0 :     watson = 0.0_dp
    1147           0 :     do i = 1, 29
    1148             : 
    1149             :        s1 = 0.0_dp
    1150             :        d = 1.0_dp
    1151           0 :        do j = 2, n
    1152           0 :           s1 = s1 + real ( j - 1, dp ) * d * x(j)
    1153           0 :           d = d * real ( i, dp ) / 29.0_dp
    1154             :        end do
    1155             : 
    1156             :        s2 = 0.0_dp
    1157             :        d = 1.0_dp
    1158           0 :        do j = 1, n
    1159           0 :           s2 = s2 + d * x(j)
    1160           0 :           d = d * real ( i, dp ) / 29.0_dp
    1161             :        end do
    1162             : 
    1163           0 :        watson = watson + ( s1 - s2 * s2 - 1.0_dp )**2
    1164             : 
    1165             :     end do
    1166             : 
    1167           0 :     watson = watson + x(1) * x(1) + ( x(2) - x(1) * x(1) - 1.0_dp )**2
    1168             : 
    1169           0 :   end function watson
    1170             : 
    1171             :   ! ------------------------------------------------------------------
    1172             :   !
    1173             :   ! The penalty function #1, 1 <= N.
    1174             :   ! Solution:
    1175             :   !    with Mathematica (numerical results)
    1176             :   !    x(1:1) = 0.5000049998750047_dp
    1177             :   !    f(x)   = 0.000002499975000499985_dp
    1178             :   !
    1179             :   !    x(1:2) = 0.35355985481744073_dp
    1180             :   !    f(x)   = 0.000008357780799989139_dp
    1181             :   !
    1182             :   !    x(1:3) = 0.2886234818387535_dp
    1183             :   !    f(x)   = 0.000015179340383244187_dp
    1184             :   !
    1185             :   !    x(1:4) = 0.2500074995875379_dp
    1186             :   !    f(x)   = 0.000022499775008999372_dp
    1187             :   !
    1188             :   !    x(1:5) = 0.2236145612000511_dp
    1189             :   !    f(x)   = 0.000030139018845277502_dp
    1190             :   !
    1191             :   !    x(1:6) = 0.20413210344548943_dp
    1192             :   !    f(x)   = 0.00003800472253975947_dp
    1193             :   !
    1194             :   !    x(1:7) = 0.18899034607915027_dp
    1195             :   !    f(x)   = 0.00004604202648884626_dp
    1196             :   !
    1197             :   !    x(1:8) = 0.1767849268724064_dp
    1198             :   !    f(x)   = 0.00005421518662591646_dp
    1199             :   !
    1200             :   !  Licensing:
    1201             :   !
    1202             :   !    This code is distributed under the GNU LGPL license.
    1203             :   !
    1204             :   !  Modified:
    1205             :   !
    1206             :   !    15 March 2000
    1207             :   !
    1208             :   !  Author:
    1209             :   !
    1210             :   !    John Burkardt
    1211             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1212             :   !
    1213             :   !  Parameters:
    1214             :   !
    1215             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1216             :   !
    1217             : 
    1218           0 :   function penalty1(x)
    1219             : 
    1220             :     implicit none
    1221             : 
    1222             :     real(dp), dimension(:), intent(in) :: x
    1223             :     real(dp)                           :: penalty1
    1224             : 
    1225             :     integer(i4)         :: n
    1226             :     real(dp), parameter :: ap = 0.00001_dp
    1227           0 :     real(dp)            :: t1
    1228           0 :     real(dp)            :: t2
    1229             : 
    1230           0 :     n = size(x)
    1231           0 :     t1 = - 0.25_dp + sum ( x(1:n)**2 )
    1232             : 
    1233           0 :     t2 = sum ( ( x(1:n) - 1.0_dp )**2 )
    1234             : 
    1235           0 :     penalty1 = ap * t2 + t1 * t1
    1236             : 
    1237           0 :   end function penalty1
    1238             : 
    1239             :   ! ------------------------------------------------------------------
    1240             :   !
    1241             :   ! The penalty function #2, 1 <= N.
    1242             :   ! Solution:
    1243             :   !    with Mathematica (numerical results)
    1244             :   !
    1245             :   !    x(1:1) = (/ 0.7914254791661984_dp /)
    1246             :   !    f(x)   = 0.4893952147007766_dp
    1247             :   !
    1248             :   !    x(1:2) = (/ 0.2000002053277478_dp, 0.95916622198851_dp /)
    1249             :   !    f(x)   = 0.000000806639004111886_dp
    1250             :   !
    1251             :   !    x(1:3) = (/ 0.19999990827773015_dp, 0.521451491987707_dp, 0.5798077888356243_dp /)
    1252             :   !    f(x)   = 0.0000031981885430064677_dp
    1253             :   !
    1254             :   !    x(1:4) = (/ 0.19999930547325262_dp, 0.19165582942317613_dp, 0.47960388408453586_dp, &
    1255             :   !                0.519390092116238_dp /)
    1256             :   !    f(x)   = 0.000009376294404291533_dp
    1257             :   !
    1258             :   !    x(1:5) = (/ 0.19999832542354226_dp, 0.09277573718478778_dp, 0.20939661885734498_dp, &
    1259             :   !                0.4467231497703405_dp, 0.4846761802898935_dp /)
    1260             :   !    f(x)   = 0.000021387590981336432_dp
    1261             :   !
    1262             :   !    x(1:6) = (/ 0.19999687534275826_dp, 0.05202606784766213_dp, 0.11181980463664613_dp, &
    1263             :   !                0.21122219156860741_dp, 0.4237150785565095_dp, 0.45116217266910164_dp /)
    1264             :   !    f(x)   = 0.00004193126194907272_dp
    1265             :   !
    1266             :   !    x(1:7) = (/ 0.1999947753716143_dp,  0.03223698574640245_dp, 0.06616491022124847_dp, &
    1267             :   !                0.11800087933141483_dp, 0.2086993295142654_dp, 0.4018822793193398_dp, &
    1268             :   !                0.4272124489678876_dp /)
    1269             :   !    f(x)   = 0.00007445577533975632_dp
    1270             :   !
    1271             :   !    x(1:8) = (/ 0.19999196971962244_dp, 0.02124196902699145_dp, 0.04164127214921979_dp, &
    1272             :   !                0.0721284757968418_dp, 0.11998769382443852_dp, 0.20359198648845672_dp, &
    1273             :   !                0.3808503545095562_dp,  0.41039238853670246_dp /)
    1274             :   !    f(x)   = 0.0001233351431976925_dp
    1275             :   !
    1276             :   !  Licensing:
    1277             :   !
    1278             :   !    This code is distributed under the GNU LGPL license.
    1279             :   !
    1280             :   !  Modified:
    1281             :   !
    1282             :   !    15 March 2000
    1283             :   !
    1284             :   !  Author:
    1285             :   !
    1286             :   !    John Burkardt
    1287             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1288             :   !
    1289             :   !  Parameters:
    1290             :   !
    1291             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1292             :   !
    1293             : 
    1294           0 :   function penalty2(x)
    1295             : 
    1296             :     implicit none
    1297             : 
    1298             :     integer(i4) :: n
    1299             : 
    1300             :     real(dp), parameter :: ap = 0.00001_dp
    1301           0 :     real(dp) :: d2
    1302             :     real(dp) :: penalty2
    1303             :     integer(i4) ::j
    1304           0 :     real(dp) :: s1
    1305           0 :     real(dp) :: s2
    1306           0 :     real(dp) :: s3
    1307           0 :     real(dp) :: t1
    1308           0 :     real(dp) :: t2
    1309           0 :     real(dp) :: t3
    1310             :     real(dp), dimension(:), intent(in) :: x
    1311             : 
    1312           0 :     n = size(x)
    1313           0 :     t1 = -1.0_dp
    1314           0 :     t2 = 0.0_dp
    1315           0 :     t3 = 0.0_dp
    1316           0 :     d2 = 1.0_dp
    1317           0 :     s2 = 0.0_dp
    1318           0 :     do j = 1, n
    1319           0 :        t1 = t1 + real ( n - j + 1, dp ) * x(j)**2
    1320           0 :        s1 = exp ( x(j) / 10.0_dp )
    1321           0 :        if ( 1 < j ) then
    1322           0 :           s3 = s1 + s2 - d2 * ( exp ( 0.1_dp ) + 1.0_dp )
    1323           0 :           t2 = t2 + s3 * s3
    1324           0 :           t3 = t3 + ( s1 - 1.0_dp / exp ( 0.1_dp ) )**2
    1325             :        end if
    1326           0 :        s2 = s1
    1327           0 :        d2 = d2 * exp ( 0.1_dp )
    1328             :     end do
    1329             : 
    1330           0 :     penalty2 = ap * ( t2 + t3 ) + t1 * t1 + ( x(1) - 0.2_dp )**2
    1331             : 
    1332           0 :   end function penalty2
    1333             : 
    1334             :   ! ------------------------------------------------------------------
    1335             :   !
    1336             :   ! The Brown badly scaled function, N = 2.
    1337             :   ! Solution: x(1:2) = (/ 1.0D+06, 2.0D-06 /)
    1338             :   !
    1339             :   !  Licensing:
    1340             :   !
    1341             :   !    This code is distributed under the GNU LGPL license.
    1342             :   !
    1343             :   !  Modified:
    1344             :   !
    1345             :   !    15 March 2000
    1346             :   !
    1347             :   !  Author:
    1348             :   !
    1349             :   !    John Burkardt
    1350             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1351             :   !
    1352             :   !  Parameters:
    1353             :   !
    1354             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1355             :   !
    1356             : 
    1357           0 :   function brown_badly_scaled(x)
    1358             : 
    1359             :     implicit none
    1360             : 
    1361             :     !    integer(i4) :: n
    1362             : 
    1363             :     real(dp) :: brown_badly_scaled
    1364             :     real(dp), dimension(:), intent(in) :: x
    1365             : 
    1366           0 :     brown_badly_scaled = ( x(1) - 1000000.0_dp )**2 &
    1367           0 :          + ( x(2) - 0.000002_dp )**2 &
    1368           0 :          + ( x(1) * x(2) - 2.0_dp )**2
    1369             : 
    1370           0 :   end function brown_badly_scaled
    1371             : 
    1372             :   ! ------------------------------------------------------------------
    1373             :   !
    1374             :   ! The Brown and Dennis function, N = 4.
    1375             :   ! Solution:
    1376             :   !    x(1:4) = (/ -11.5944399230_dp, 13.2036300657_dp, &
    1377             :   !                -0.4034395329_dp,   0.2367787297_dp /)
    1378             :   !
    1379             :   !  Licensing:
    1380             :   !
    1381             :   !    This code is distributed under the GNU LGPL license.
    1382             :   !
    1383             :   !  Modified:
    1384             :   !
    1385             :   !    05 May 2000
    1386             :   !
    1387             :   !  Author:
    1388             :   !
    1389             :   !    John Burkardt
    1390             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1391             :   !
    1392             :   !  Parameters:
    1393             :   !
    1394             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1395             :   !
    1396             : 
    1397           0 :   function brown_dennis(x)
    1398             : 
    1399             :     implicit none
    1400             : 
    1401             :     !    integer(i4) :: n
    1402             : 
    1403           0 :     real(dp) :: c
    1404             :     real(dp) :: brown_dennis
    1405           0 :     real(dp) :: f1
    1406           0 :     real(dp) :: f2
    1407             :     integer(i4) ::i
    1408             :     real(dp), dimension(:), intent(in) :: x
    1409             : 
    1410           0 :     brown_dennis = 0.0_dp
    1411             : 
    1412           0 :     do i = 1, 20
    1413             : 
    1414           0 :        c = real ( i, dp ) / 5.0_dp
    1415           0 :        f1 = x(1) + c * x(2) - exp ( c )
    1416           0 :        f2 = x(3) + sin ( c ) * x(4) - cos ( c )
    1417             : 
    1418           0 :        brown_dennis = brown_dennis + f1**4 + 2.0_dp * f1 * f1 * f2 * f2 + f2**4
    1419             : 
    1420             :     end do
    1421             : 
    1422           0 :     brown_dennis = brown_dennis - 85822.2016263563_dp
    1423             : 
    1424           0 :   end function brown_dennis
    1425             : 
    1426             :   ! ------------------------------------------------------------------
    1427             :   !
    1428             :   ! The Gulf R&D function, N = 3.
    1429             :   ! Solution: x(1:3) = (/ 50.0_dp, 25.0_dp, 1.5_dp /)
    1430             :   !
    1431             :   !  Licensing:
    1432             :   !
    1433             :   !    This code is distributed under the GNU LGPL license.
    1434             :   !
    1435             :   !  Modified:
    1436             :   !
    1437             :   !    15 March 2000
    1438             :   !
    1439             :   !  Author:
    1440             :   !
    1441             :   !    John Burkardt
    1442             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1443             :   !
    1444             :   !  Parameters:
    1445             :   !
    1446             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1447             :   !
    1448             : 
    1449           0 :   function gulf_rd(x)
    1450             : 
    1451             :     implicit none
    1452             : 
    1453             :     !    integer(i4) :: n
    1454             : 
    1455           0 :     real(dp) :: arg
    1456             :     real(dp) :: gulf_rd
    1457             :     integer(i4) ::i
    1458           0 :     real(dp) :: r
    1459           0 :     real(dp) :: t
    1460             :     real(dp), dimension(:), intent(in) :: x
    1461             :     real(dp) :: sqrtHuge
    1462             : 
    1463           0 :     sqrtHuge = Huge(1.0_dp)**0.5_dp -1.0_dp
    1464             : 
    1465           0 :     gulf_rd = 0.0_dp
    1466           0 :     do i = 1, 99
    1467             : 
    1468           0 :        arg = real ( i, dp ) / 100.0_dp
    1469             :        r = abs ( ( - 50.0_dp * log ( arg ) )**( 2.0_dp / 3.0_dp ) &
    1470           0 :             + 25.0_dp - x(2) )
    1471             : 
    1472             :        ! print*, 'arg         = ',arg
    1473             :        ! print*, 'r           = ',r
    1474             :        ! print*, 'x(1)        = ',x(1)
    1475             :        ! print*, 'x(2)        = ',x(2)
    1476             :        ! print*, 'x(3)        = ',x(3)
    1477             : 
    1478             :        ! avoiding underflow
    1479           0 :        if ( x(3)*Log(r) .gt. -708.-dp ) then
    1480           0 :           print*, '-exp(x(3)*Log(r)) = ',-exp(x(3)*Log(r))
    1481           0 :           print*, 'x(1)              = ',x(1)
    1482           0 :           if ( -exp(x(3)*Log(r)) / x(1) .lt. -708._dp) then
    1483           0 :              t = -arg
    1484             :           else
    1485           0 :              if ( -r**x(3) / x(1) .gt. 708._dp) then
    1486           0 :                 t = 1000000._dp - arg
    1487             :              else
    1488           0 :                 t = exp ( - r**x(3) / x(1) ) - arg
    1489             :              end if
    1490             :           end if
    1491             :        else
    1492           0 :           t = -arg
    1493             :        end if
    1494             : 
    1495           0 :        if ( abs(t) .gt. sqrtHuge ) then
    1496             :           ! under/overflow case
    1497           0 :           t = sqrtHuge
    1498             :        end if
    1499             : 
    1500           0 :        if ( Huge(1.0_dp) -gulf_rd .gt. t*t ) then
    1501             :           ! usual case
    1502           0 :           gulf_rd = gulf_rd + t * t
    1503             :        else
    1504             :           ! overflow case
    1505             :           gulf_rd = Huge(1.0_dp)
    1506             :        end if
    1507             : 
    1508             :     end do
    1509             : 
    1510           0 :   end function gulf_rd
    1511             : 
    1512             :   ! ------------------------------------------------------------------
    1513             :   !
    1514             :   ! The trigonometric function, 1 <= N.
    1515             :   ! Solution: x(1:n) = 0.0_dp
    1516             :   !
    1517             :   !  Licensing:
    1518             :   !
    1519             :   !    This code is distributed under the GNU LGPL license.
    1520             :   !
    1521             :   !  Modified:
    1522             :   !
    1523             :   !    15 March 2000
    1524             :   !
    1525             :   !  Author:
    1526             :   !
    1527             :   !    John Burkardt
    1528             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1529             :   !
    1530             :   !  Parameters:
    1531             :   !
    1532             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1533             :   !
    1534             : 
    1535           0 :   function trigonometric(x)
    1536             : 
    1537             :     implicit none
    1538             : 
    1539             :     integer(i4) :: n
    1540             : 
    1541             :     real(dp) :: trigonometric
    1542             :     integer(i4) ::j
    1543           0 :     real(dp) :: s1
    1544           0 :     real(dp) :: t
    1545             :     real(dp), dimension(:), intent(in) :: x
    1546             : 
    1547           0 :     n = size(x)
    1548           0 :     s1 = sum ( cos ( x(1:n) ) )
    1549             : 
    1550             :     trigonometric = 0.0_dp
    1551           0 :     do j = 1, n
    1552           0 :        t = real ( n + j, dp ) - sin ( x(j) ) &
    1553           0 :             - s1 - real ( j, dp ) * cos ( x(j) )
    1554           0 :        trigonometric = trigonometric + t * t
    1555             :     end do
    1556             : 
    1557           0 :   end function trigonometric
    1558             : 
    1559             :   ! ------------------------------------------------------------------
    1560             :   !
    1561             :   ! The extended Rosenbrock parabolic valley function, 1 <= N.
    1562             :   ! Solution: x(1:n) = 1.0_dp
    1563             :   !
    1564             :   !  Licensing:
    1565             :   !
    1566             :   !    This code is distributed under the GNU LGPL license.
    1567             :   !
    1568             :   !  Modified:
    1569             :   !
    1570             :   !    15 March 2000
    1571             :   !
    1572             :   !  Author:
    1573             :   !
    1574             :   !    John Burkardt
    1575             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1576             :   !
    1577             :   !  Reference:
    1578             :   !
    1579             :   !    Richard Brent,
    1580             :   !    Algorithms for Minimization with Derivatives,
    1581             :   !    Dover, 2002,
    1582             :   !    ISBN: 0-486-41998-3,
    1583             :   !    LC: QA402.5.B74.
    1584             :   !
    1585             :   !  Parameters:
    1586             :   !
    1587             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1588             :   !
    1589             : 
    1590           0 :   function ext_rosenbrock_parabolic_valley(x)
    1591             : 
    1592             :     implicit none
    1593             : 
    1594             :     integer(i4) :: n
    1595             : 
    1596             :     real(dp) :: ext_rosenbrock_parabolic_valley
    1597             :     integer(i4) ::j
    1598             :     real(dp), dimension(:), intent(in) :: x
    1599             : 
    1600           0 :     n = size(x)
    1601           0 :     ext_rosenbrock_parabolic_valley = 0.0_dp
    1602           0 :     do j = 1, n
    1603           0 :        if ( mod ( j, 2 ) == 1 ) then
    1604           0 :           ext_rosenbrock_parabolic_valley = ext_rosenbrock_parabolic_valley + ( 1.0_dp - x(j) )**2
    1605             :        else
    1606           0 :           ext_rosenbrock_parabolic_valley = ext_rosenbrock_parabolic_valley + 100.0_dp * ( x(j) - x(j-1) * x(j-1) )**2
    1607             :        end if
    1608             :     end do
    1609             : 
    1610           0 :   end function ext_rosenbrock_parabolic_valley
    1611             : 
    1612             :   ! ------------------------------------------------------------------
    1613             :   !
    1614             :   ! The extended Powell singular quartic function, 4 <= N.
    1615             :   ! Solution: x(1:n) = 0.0_dp
    1616             :   !
    1617             :   !  Discussion:
    1618             :   !
    1619             :   !    The Hessian matrix is doubly singular at the minimizer,
    1620             :   !    suggesting that most optimization routines will experience
    1621             :   !    a severe slowdown in convergence.
    1622             :   !
    1623             :   !    The problem is usually only defined for N being a multiple of 4.
    1624             :   !
    1625             :   !  Licensing:
    1626             :   !
    1627             :   !    This code is distributed under the GNU LGPL license.
    1628             :   !
    1629             :   !  Modified:
    1630             :   !
    1631             :   !    05 May 2000
    1632             :   !
    1633             :   !  Author:
    1634             :   !
    1635             :   !    John Burkardt
    1636             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1637             :   !
    1638             :   !  Reference:
    1639             :   !
    1640             :   !    Richard Brent,
    1641             :   !    Algorithms for Minimization with Derivatives,
    1642             :   !    Dover, 2002,
    1643             :   !    ISBN: 0-486-41998-3,
    1644             :   !    LC: QA402.5.B74.
    1645             :   !
    1646             :   !  Parameters:
    1647             :   !
    1648             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1649             :   !
    1650             : 
    1651           0 :   function ext_powell_singular_quartic(x)
    1652             : 
    1653             :     implicit none
    1654             : 
    1655             :     integer(i4) :: n
    1656             : 
    1657             :     real(dp) :: ext_powell_singular_quartic
    1658           0 :     real(dp) :: f1
    1659           0 :     real(dp) :: f2
    1660           0 :     real(dp) :: f3
    1661           0 :     real(dp) :: f4
    1662             :     integer(i4) ::j
    1663             :     real(dp), dimension(:), intent(in) :: x
    1664           0 :     real(dp) :: xjp1
    1665           0 :     real(dp) :: xjp2
    1666           0 :     real(dp) :: xjp3
    1667             : 
    1668           0 :     n = size(x)
    1669           0 :     ext_powell_singular_quartic = 0.0_dp
    1670             : 
    1671           0 :     do j = 1, n, 4
    1672             : 
    1673           0 :        if ( j + 1 <= n ) then
    1674           0 :           xjp1 = x(j+1)
    1675             :        else
    1676             :           xjp1 = 0.0_dp
    1677             :        end if
    1678             : 
    1679           0 :        if ( j + 2 <= n ) then
    1680           0 :           xjp2 = x(j+2)
    1681             :        else
    1682             :           xjp2 = 0.0_dp
    1683             :        end if
    1684             : 
    1685           0 :        if ( j + 3 <= n ) then
    1686           0 :           xjp3 = x(j+3)
    1687             :        else
    1688             :           xjp3 = 0.0_dp
    1689             :        end if
    1690             : 
    1691           0 :        f1 = x(j) + 10.0_dp * xjp1
    1692             : 
    1693           0 :        if ( j + 1 <= n ) then
    1694           0 :           f2 = xjp2 - xjp3
    1695             :        else
    1696             :           f2 = 0.0_dp
    1697             :        end if
    1698             : 
    1699           0 :        if ( j + 2 <= n ) then
    1700           0 :           f3 = xjp1 - 2.0_dp * xjp2
    1701             :        else
    1702             :           f3 = 0.0_dp
    1703             :        end if
    1704             : 
    1705           0 :        if ( j + 3 <= n ) then
    1706           0 :           f4 = x(j) - xjp3
    1707             :        else
    1708             :           f4 = 0.0_dp
    1709             :        end if
    1710             : 
    1711             :        ext_powell_singular_quartic = ext_powell_singular_quartic +            f1 * f1 &
    1712             :             +  5.0_dp * f2 * f2 &
    1713             :             +            f3 * f3 * f3 * f3 &
    1714           0 :             + 10.0_dp * f4 * f4 * f4 * f4
    1715             : 
    1716             :     end do
    1717             : 
    1718           0 :   end function ext_powell_singular_quartic
    1719             : 
    1720             :   ! ------------------------------------------------------------------
    1721             :   !
    1722             :   ! The Beale function, N = 2.
    1723             :   ! Solution: x(1:2) = (/ 3.0_dp, 0.5_dp /)
    1724             :   ! Search domain: -4.5 <= xi <= 4.5, i = 1, 2.
    1725             :   !
    1726             :   !  Discussion:
    1727             :   !
    1728             :   !    Range according to:
    1729             :   !       http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/
    1730             :   !       Hedar_files/TestGO_files/Page288.htm
    1731             :   !
    1732             :   !    This function has a valley approaching the line X(2) = 1.
    1733             :   !
    1734             :   !    The function has a global minimizer:
    1735             :   !
    1736             :   !      X(*) = ( 3.0, 0.5 ), F(X*) = 0.0.
    1737             :   !
    1738             :   !  Licensing:
    1739             :   !
    1740             :   !    This code is distributed under the GNU LGPL license.
    1741             :   !
    1742             :   !  Modified:
    1743             :   !
    1744             :   !    28 January 2008
    1745             :   !
    1746             :   !  Author:
    1747             :   !
    1748             :   !    John Burkardt
    1749             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1750             :   !
    1751             :   !  Reference:
    1752             :   !
    1753             :   !    Evelyn Beale,
    1754             :   !    On an Iterative Method for Finding a Local Minimum of a Function
    1755             :   !    of More than One Variable,
    1756             :   !    Technical Report 25,
    1757             :   !    Statistical Techniques Research Group,
    1758             :   !    Princeton University, 1958.
    1759             :   !
    1760             :   !    Richard Brent,
    1761             :   !    Algorithms for Minimization with Derivatives,
    1762             :   !    Dover, 2002,
    1763             :   !    ISBN: 0-486-41998-3,
    1764             :   !    LC: QA402.5.B74.
    1765             :   !
    1766             :   !  Parameters:
    1767             :   !
    1768             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1769             :   !
    1770             : 
    1771           0 :   function beale(x)
    1772             : 
    1773             :     implicit none
    1774             : 
    1775             :     !    integer(i4) :: n
    1776             : 
    1777             :     real(dp) :: beale
    1778           0 :     real(dp) :: f1
    1779           0 :     real(dp) :: f2
    1780           0 :     real(dp) :: f3
    1781             :     real(dp), dimension(:), intent(in) :: x
    1782             : 
    1783           0 :     f1 = 1.5_dp   - x(1) * ( 1.0_dp - x(2)    )
    1784           0 :     f2 = 2.25_dp  - x(1) * ( 1.0_dp - x(2) * x(2) )
    1785           0 :     f3 = 2.625_dp - x(1) * ( 1.0_dp - x(2) * x(2) * x(2) )
    1786             : 
    1787           0 :     beale = f1 * f1 + f2 * f2 + f3 * f3
    1788             : 
    1789           0 :   end function beale
    1790             : 
    1791             :   ! ------------------------------------------------------------------
    1792             :   !
    1793             :   ! The Wood function, N = 4.
    1794             :   ! Solution: x(1:4) = (/ 1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp /)
    1795             :   !
    1796             :   !  Licensing:
    1797             :   !
    1798             :   !    This code is distributed under the GNU LGPL license.
    1799             :   !
    1800             :   !  Modified:
    1801             :   !
    1802             :   !    06 January 2008
    1803             :   !
    1804             :   !  Author:
    1805             :   !
    1806             :   !    John Burkardt
    1807             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1808             :   !
    1809             :   !  Reference:
    1810             :   !
    1811             :   !    Richard Brent,
    1812             :   !    Algorithms for Minimization with Derivatives,
    1813             :   !    Dover, 2002,
    1814             :   !    ISBN: 0-486-41998-3,
    1815             :   !    LC: QA402.5.B74.
    1816             :   !
    1817             :   !  Parameters:
    1818             :   !
    1819             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1820             :   !
    1821             : 
    1822           0 :   function wood(x)
    1823             : 
    1824             :     implicit none
    1825             : 
    1826             :     !    integer(i4) :: n
    1827             : 
    1828             :     real(dp) :: wood
    1829           0 :     real(dp) :: f1
    1830           0 :     real(dp) :: f2
    1831           0 :     real(dp) :: f3
    1832           0 :     real(dp) :: f4
    1833           0 :     real(dp) :: f5
    1834           0 :     real(dp) :: f6
    1835             :     real(dp), dimension(:), intent(in) :: x
    1836             : 
    1837           0 :     f1 = x(2) - x(1) * x(1)
    1838           0 :     f2 = 1.0_dp - x(1)
    1839           0 :     f3 = x(4) - x(3) * x(3)
    1840           0 :     f4 = 1.0_dp - x(3)
    1841           0 :     f5 = x(2) + x(4) - 2.0_dp
    1842           0 :     f6 = x(2) - x(4)
    1843             : 
    1844             :     wood = 100.0_dp * f1 * f1 &
    1845             :          +             f2 * f2 &
    1846             :          +  90.0_dp * f3 * f3 &
    1847             :          +             f4 * f4 &
    1848             :          +  10.0_dp * f5 * f5 &
    1849           0 :          +   0.1_dp * f6 * f6
    1850             : 
    1851           0 :   end function wood
    1852             : 
    1853             :   ! ------------------------------------------------------------------
    1854             :   !
    1855             :   ! The Chebyquad function, 1 <= N.
    1856             :   ! Solution: n==2: x(1:2) = (/ 0.2113249_dp, 0.7886751_dp /)
    1857             :   !           n==4: x(1:4) = (/ 0.1026728_dp, 0.4062037_dp, 0.5937963_dp, 0.8973272_dp /)
    1858             :   !           n==6: x(1:6) = (/ 0.066877_dp, 0.288741_dp, 0.366682_dp,
    1859             :   !                             0.633318_dp, 0.711259_dp, 0.933123_dp /)
    1860             :   !           n==8: x(1:8) = (/ 0.043153_dp, 0.193091_dp, 0.266329_dp, 0.500000_dp, &
    1861             :   !                             0.500000_dp, 0.733671_dp, 0.806910_dp, 0.956847_dp /)
    1862             :   !           else unknown
    1863             :   !
    1864             :   !  Licensing:
    1865             :   !
    1866             :   !    This code is distributed under the GNU LGPL license.
    1867             :   !
    1868             :   !  Modified:
    1869             :   !
    1870             :   !    23 March 2000
    1871             :   !
    1872             :   !  Author:
    1873             :   !
    1874             :   !    John Burkardt
    1875             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1876             :   !
    1877             :   !  Reference:
    1878             :   !
    1879             :   !    Richard Brent,
    1880             :   !    Algorithms for Minimization with Derivatives,
    1881             :   !    Dover, 2002,
    1882             :   !    ISBN: 0-486-41998-3,
    1883             :   !    LC: QA402.5.B74.
    1884             :   !
    1885             :   !  Parameters:
    1886             :   !
    1887             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1888             :   !
    1889             : 
    1890           0 :   function chebyquad(x)
    1891             : 
    1892             :     implicit none
    1893             : 
    1894             :     integer(i4) :: n
    1895             : 
    1896             :     real(dp) :: chebyquad
    1897             :     real(dp), dimension(:), intent(in) :: x
    1898             : 
    1899           0 :     real(dp), dimension(size(x)) :: fvec
    1900             :     integer(i4) :: i
    1901             :     integer(i4) :: j
    1902           0 :     real(dp) :: t
    1903           0 :     real(dp) :: t1
    1904           0 :     real(dp) :: t2
    1905           0 :     real(dp) :: th
    1906             : 
    1907             :     !
    1908             :     !  Compute FVEC.
    1909             :     !
    1910           0 :     n = size(x)
    1911           0 :     fvec(1:n) = 0.0_dp
    1912           0 :     do j = 1, n
    1913           0 :        t1 = 1.0_dp
    1914           0 :        t2 = 2.0_dp * x(j) - 1.0_dp
    1915           0 :        t = 2.0_dp * t2
    1916           0 :        do i = 1, n
    1917           0 :           fvec(i) = fvec(i) + t2
    1918           0 :           th = t * t2 - t1
    1919           0 :           t1 = t2
    1920           0 :           t2 = th
    1921             :        end do
    1922             :     end do
    1923             : 
    1924           0 :     do i = 1, n
    1925           0 :        fvec(i) = fvec(i) / real ( n, dp )
    1926           0 :        if ( mod ( i, 2 ) == 0 ) then
    1927           0 :           fvec(i) = fvec(i) + 1.0_dp / ( real ( i, dp )**2 - 1.0_dp )
    1928             :        end if
    1929             :     end do
    1930             :     !call p18_fvec ( n, x, fvec )
    1931             :     !
    1932             :     !  Compute F.
    1933             :     !
    1934           0 :     chebyquad = sum ( fvec(1:n)**2 )
    1935             : 
    1936           0 :   end function chebyquad
    1937             : 
    1938             :   ! ------------------------------------------------------------------
    1939             :   !
    1940             :   ! The Leon''s cubic valley function, N = 2.
    1941             :   ! Solution: x(1:2) = (/ 1.0_dp, 1.0_dp /)
    1942             :   !
    1943             :   !  Discussion:
    1944             :   !
    1945             :   !    The function is similar to Rosenbrock's.  The "valley" follows
    1946             :   !    the curve X(2) = X(1)**3.
    1947             :   !
    1948             :   !  Licensing:
    1949             :   !
    1950             :   !    This code is distributed under the GNU LGPL license.
    1951             :   !
    1952             :   !  Modified:
    1953             :   !
    1954             :   !    17 March 2000
    1955             :   !
    1956             :   !  Author:
    1957             :   !
    1958             :   !    John Burkardt
    1959             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    1960             :   !
    1961             :   !  Reference:
    1962             :   !
    1963             :   !    Richard Brent,
    1964             :   !    Algorithms for Minimization with Derivatives,
    1965             :   !    Dover, 2002,
    1966             :   !    ISBN: 0-486-41998-3,
    1967             :   !    LC: QA402.5.B74.
    1968             :   !
    1969             :   !    A Leon,
    1970             :   !    A Comparison of Eight Known Optimizing Procedures,
    1971             :   !    in Recent Advances in Optimization Techniques,
    1972             :   !    edited by Abraham Lavi, Thomas Vogl,
    1973             :   !    Wiley, 1966.
    1974             :   !
    1975             :   !  Parameters:
    1976             :   !
    1977             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    1978             :   !
    1979             : 
    1980           0 :   function leon_cubic_valley(x)
    1981             : 
    1982             :     implicit none
    1983             : 
    1984             :     !    integer(i4) :: n
    1985             : 
    1986             :     real(dp) :: leon_cubic_valley
    1987           0 :     real(dp) :: f1
    1988           0 :     real(dp) :: f2
    1989             :     real(dp), dimension(:), intent(in) :: x
    1990             : 
    1991           0 :     f1 = x(2) - x(1) * x(1) * x(1)
    1992           0 :     f2 = 1.0_dp - x(1)
    1993             : 
    1994             :     leon_cubic_valley = 100.0_dp * f1 * f1 &
    1995           0 :          +             f2 * f2
    1996             : 
    1997           0 :   end function leon_cubic_valley
    1998             : 
    1999             :   ! ------------------------------------------------------------------
    2000             :   !
    2001             :   ! The Gregory and Karney''s Tridiagonal Matrix Function, 1 <= N.
    2002             :   ! Solution: forall(i=1:n) x(i) = real(n+1-i, dp)
    2003             :   !
    2004             :   !  Discussion:
    2005             :   !
    2006             :   !    The function has the form
    2007             :   !      f = x'*A*x - 2*x(1)
    2008             :   !    where A is the (-1,2,-1) tridiagonal matrix, except that A(1,1) is 1.
    2009             :   !    The minimum value of F(X) is -N, which occurs for
    2010             :   !      x = ( n, n-1, ..., 2, 1 ).
    2011             :   !
    2012             :   !  Licensing:
    2013             :   !
    2014             :   !    This code is distributed under the GNU LGPL license.
    2015             :   !
    2016             :   !  Modified:
    2017             :   !
    2018             :   !    20 March 2000
    2019             :   !
    2020             :   !  Author:
    2021             :   !
    2022             :   !    John Burkardt
    2023             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2024             :   !
    2025             :   !  Reference:
    2026             :   !
    2027             :   !    Richard Brent,
    2028             :   !    Algorithms for Minimization with Derivatives,
    2029             :   !    Prentice Hall, 1973,
    2030             :   !    Reprinted by Dover, 2002.
    2031             :   !
    2032             :   !  Parameters:
    2033             :   !
    2034             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2035             :   !
    2036             : 
    2037           0 :   function gregory_karney_tridia_matrix(x)
    2038             : 
    2039             :     implicit none
    2040             : 
    2041             :     integer(i4) :: n
    2042             : 
    2043             :     real(dp) :: gregory_karney_tridia_matrix
    2044             :     integer(i4) ::i
    2045             :     real(dp), dimension(:), intent(in) :: x
    2046             : 
    2047           0 :     n = size(x)
    2048           0 :     gregory_karney_tridia_matrix = x(1) * x(1) + 2.0_dp * sum ( x(2:n)**2 )
    2049             : 
    2050           0 :     do i = 1, n-1
    2051           0 :        gregory_karney_tridia_matrix = gregory_karney_tridia_matrix - 2.0_dp * x(i) * x(i+1)
    2052             :     end do
    2053             : 
    2054           0 :     gregory_karney_tridia_matrix = gregory_karney_tridia_matrix - 2.0_dp * x(1)
    2055             : 
    2056           0 :     gregory_karney_tridia_matrix = gregory_karney_tridia_matrix + real(n,dp)
    2057             : 
    2058           0 :   end function gregory_karney_tridia_matrix
    2059             : 
    2060             :   ! ------------------------------------------------------------------
    2061             :   !
    2062             :   ! The Hilbert function, 1 <= N.
    2063             :   ! Solution: x(1:n) = 0.0_dp
    2064             :   !
    2065             :   !  Discussion:
    2066             :   !
    2067             :   !    The function has the form
    2068             :   !      f = x'*A*x
    2069             :   !    where A is the Hilbert matrix.  The minimum value
    2070             :   !    of F(X) is 0, which occurs for
    2071             :   !      x = ( 0, 0, ..., 0 ).
    2072             :   !
    2073             :   !  Licensing:
    2074             :   !
    2075             :   !    This code is distributed under the GNU LGPL license.
    2076             :   !
    2077             :   !  Modified:
    2078             :   !
    2079             :   !    20 March 2000
    2080             :   !
    2081             :   !  Author:
    2082             :   !
    2083             :   !    John Burkardt
    2084             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2085             :   !
    2086             :   !  Reference:
    2087             :   !
    2088             :   !    Richard Brent,
    2089             :   !    Algorithms for Minimization with Derivatives,
    2090             :   !    Dover, 2002,
    2091             :   !    ISBN: 0-486-41998-3,
    2092             :   !    LC: QA402.5.B74.
    2093             :   !
    2094             :   !  Parameters:
    2095             :   !
    2096             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2097             :   !
    2098             : 
    2099           0 :   function hilbert(x)
    2100             : 
    2101             :     implicit none
    2102             : 
    2103             :     integer(i4) :: n
    2104             : 
    2105             :     real(dp) :: hilbert
    2106             :     integer(i4) ::i
    2107             :     integer(i4) ::j
    2108             :     real(dp), dimension(:), intent(in) :: x
    2109             : 
    2110           0 :     n = size(x)
    2111           0 :     hilbert = 0.0_dp
    2112             : 
    2113           0 :     do i = 1, n
    2114           0 :        do j = 1, n
    2115           0 :           hilbert = hilbert + x(i) * x(j) / real ( i + j - 1, dp )
    2116             :        end do
    2117             :     end do
    2118             : 
    2119           0 :   end function hilbert
    2120             : 
    2121             :   ! ------------------------------------------------------------------
    2122             :   !
    2123             :   ! The De Jong Function F1, N = 3.
    2124             :   ! Solution: x(1:n) = 0.0_dp
    2125             :   !
    2126             :   !  Licensing:
    2127             :   !
    2128             :   !    This code is distributed under the GNU LGPL license.
    2129             :   !
    2130             :   !  Modified:
    2131             :   !
    2132             :   !    30 December 2000
    2133             :   !
    2134             :   !  Author:
    2135             :   !
    2136             :   !    John Burkardt
    2137             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2138             :   !
    2139             :   !  Reference:
    2140             :   !
    2141             :   !    Zbigniew Michalewicz,
    2142             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2143             :   !    Third Edition,
    2144             :   !    Springer Verlag, 1996,
    2145             :   !    ISBN: 3-540-60676-9,
    2146             :   !    LC: QA76.618.M53.
    2147             :   !
    2148             :   !  Parameters:
    2149             :   !
    2150             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2151             :   !
    2152             : 
    2153           0 :   function de_jong_f1(x)
    2154             : 
    2155             :     implicit none
    2156             : 
    2157             :     integer(i4) :: n
    2158             : 
    2159             :     real(dp) :: de_jong_f1
    2160             :     real(dp), dimension(:), intent(in) :: x
    2161             : 
    2162           0 :     n = size(x)
    2163           0 :     de_jong_f1 = sum ( x(1:n)**2 )
    2164             : 
    2165           0 :   end function de_jong_f1
    2166             : 
    2167             :   ! ------------------------------------------------------------------
    2168             :   !
    2169             :   ! The De Jong Function F2, N = 2.
    2170             :   ! Solution: x(1:n) = 1.0_dp
    2171             :   !
    2172             :   !  Licensing:
    2173             :   !
    2174             :   !    This code is distributed under the GNU LGPL license.
    2175             :   !
    2176             :   !  Modified:
    2177             :   !
    2178             :   !    31 December 2000
    2179             :   !
    2180             :   !  Author:
    2181             :   !
    2182             :   !    John Burkardt
    2183             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2184             :   !
    2185             :   !  Reference:
    2186             :   !
    2187             :   !    Zbigniew Michalewicz,
    2188             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2189             :   !    Third Edition,
    2190             :   !    Springer Verlag, 1996,
    2191             :   !    ISBN: 3-540-60676-9,
    2192             :   !    LC: QA76.618.M53.
    2193             :   !
    2194             :   !  Parameters:
    2195             :   !
    2196             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2197             :   !
    2198             : 
    2199           0 :   function de_jong_f2(x)
    2200             : 
    2201             :     implicit none
    2202             : 
    2203             :     !    integer(i4) :: n
    2204             : 
    2205             :     real(dp) :: de_jong_f2
    2206             :     real(dp), dimension(:), intent(in) :: x
    2207             : 
    2208           0 :     de_jong_f2 = 100.0_dp * ( x(1) * x(1) - x(2) )**2 + ( 1.0_dp - x(1) )**2
    2209             : 
    2210           0 :   end function de_jong_f2
    2211             : 
    2212             :   ! ------------------------------------------------------------------
    2213             :   !
    2214             :   ! The De Jong Function F3 (discontinuous), N = 5.
    2215             :   ! Solution:
    2216             :   !    x(1:n) depends on search space
    2217             :   !           = sum of integer part of left boundary values
    2218             :   !
    2219             :   !  Licensing:
    2220             :   !
    2221             :   !    This code is distributed under the GNU LGPL license.
    2222             :   !
    2223             :   !  Modified:
    2224             :   !
    2225             :   !    31 December 2000
    2226             :   !
    2227             :   !  Author:
    2228             :   !
    2229             :   !    John Burkardt
    2230             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2231             :   !
    2232             :   !  Reference:
    2233             :   !
    2234             :   !    Zbigniew Michalewicz,
    2235             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2236             :   !    Third Edition,
    2237             :   !    Springer Verlag, 1996,
    2238             :   !    ISBN: 3-540-60676-9,
    2239             :   !    LC: QA76.618.M53.
    2240             :   !
    2241             :   !  Parameters:
    2242             :   !
    2243             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2244             :   !
    2245             : 
    2246           0 :   function de_jong_f3(x)
    2247             : 
    2248             :     implicit none
    2249             : 
    2250             :     integer(i4) :: n
    2251             : 
    2252             :     real(dp) :: de_jong_f3
    2253             :     real(dp), dimension(:), intent(in) :: x
    2254             : 
    2255           0 :     n = size(x)
    2256             :     ! Original: conversion to int only possible up to i8, else "Floating invalid operation - aborting"
    2257             :     ! de_jong_f3 = real ( sum ( int ( x(1:n) ) ), dp )
    2258             : 
    2259           0 :     de_jong_f3 = sum ( aint ( x(1:n) ) )
    2260             : 
    2261           0 :   end function de_jong_f3
    2262             : 
    2263             :   ! ------------------------------------------------------------------
    2264             :   !
    2265             :   ! The De Jong Function F4 (Gaussian noise), N = 30.
    2266             :   ! Solution: x(1:n) = 0.0_dp
    2267             :   !
    2268             :   !  Discussion:
    2269             :   !
    2270             :   !    The function includes Gaussian noise, multiplied by a parameter P.
    2271             :   !
    2272             :   !    If P is zero, then the function is a proper function, and evaluating
    2273             :   !    it twice with the same argument will yield the same results.
    2274             :   !    Moreover, P25_G and P25_H are the correct gradient and hessian routines.
    2275             :   !
    2276             :   !    If P is nonzero, then evaluating the function at the same argument
    2277             :   !    twice will generally yield two distinct values; this means the function
    2278             :   !    is not even a well defined mathematical function, let alone continuous;
    2279             :   !    the gradient and hessian are not correct.  And yet, at least for small
    2280             :   !    values of P, it may be possible to approximate the minimizer of the
    2281             :   !    (underlying well-defined ) function.
    2282             :   !
    2283             :   !    The value of the parameter P is by default 1.  The user can manipulate
    2284             :   !    this value by calling P25_P_GET or P25_P_SET.
    2285             :   !
    2286             :   !  Licensing:
    2287             :   !
    2288             :   !    This code is distributed under the GNU LGPL license.
    2289             :   !
    2290             :   !  Modified:
    2291             :   !
    2292             :   !    22 January 2001
    2293             :   !
    2294             :   !  Author:
    2295             :   !
    2296             :   !    John Burkardt
    2297             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2298             :   !
    2299             :   !  Reference:
    2300             :   !
    2301             :   !    Zbigniew Michalewicz,
    2302             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2303             :   !    Third Edition,
    2304             :   !    Springer Verlag, 1996,
    2305             :   !    ISBN: 3-540-60676-9,
    2306             :   !    LC: QA76.618.M53.
    2307             :   !
    2308             :   !  Parameters:
    2309             :   !
    2310             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2311             :   !
    2312             : 
    2313           0 :   function de_jong_f4(x)
    2314             : 
    2315             :     implicit none
    2316             : 
    2317             :     integer(i4) :: n
    2318             : 
    2319             :     real(dp) :: de_jong_f4
    2320           0 :     real(dp) :: gauss
    2321             :     integer(i4) ::i
    2322           0 :     real(dp) :: p
    2323             :     real(dp), dimension(:), intent(in) :: x
    2324             : 
    2325             :     real(dp), save :: p_save = 1.0_dp
    2326             : 
    2327           0 :     n = size(x)
    2328           0 :     p = p_save
    2329             :     !call p25_p_get ( p )
    2330             : 
    2331           0 :     call normal_01_sample ( gauss )
    2332             : 
    2333           0 :     de_jong_f4 = p * gauss
    2334           0 :     do i = 1, n
    2335           0 :        de_jong_f4 = de_jong_f4 + real ( i, dp ) * x(i)**4
    2336             :     end do
    2337             : 
    2338           0 :   end function de_jong_f4
    2339             : 
    2340             :   ! ------------------------------------------------------------------
    2341             :   !
    2342             :   ! The De Jong Function F5, N = 2.
    2343             :   ! Solution: x(1:n) = 0.0_dp
    2344             :   !
    2345             :   !  Licensing:
    2346             :   !
    2347             :   !    This code is distributed under the GNU LGPL license.
    2348             :   !
    2349             :   !  Modified:
    2350             :   !
    2351             :   !    01 January 2001
    2352             :   !
    2353             :   !  Author:
    2354             :   !
    2355             :   !    John Burkardt
    2356             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2357             :   !
    2358             :   !  Reference:
    2359             :   !
    2360             :   !    Zbigniew Michalewicz,
    2361             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2362             :   !    Third Edition,
    2363             :   !    Springer Verlag, 1996,
    2364             :   !    ISBN: 3-540-60676-9,
    2365             :   !    LC: QA76.618.M53.
    2366             :   !
    2367             :   !  Parameters:
    2368             :   !
    2369             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2370             :   !
    2371             : 
    2372           0 :   function de_jong_f5(x)
    2373             : 
    2374             :     implicit none
    2375             : 
    2376             :     !    integer(i4) :: n
    2377             : 
    2378             :     integer(i4) ::a1
    2379             :     integer(i4) ::a2
    2380             :     real(dp) :: de_jong_f5
    2381             :     real(dp) :: fi
    2382           0 :     real(dp) :: fj
    2383             :     integer(i4) ::j
    2384             :     integer(i4) ::j1
    2385             :     integer(i4) ::j2
    2386             :     integer(i4), parameter :: jroot = 5
    2387             :     integer(i4), parameter :: k = 500
    2388             :     real(dp), dimension(:), intent(in) :: x
    2389             : 
    2390           0 :     fi = real(k,dp)
    2391             : 
    2392           0 :     do j=1, jroot * jroot
    2393             : 
    2394           0 :        j1 = mod(j-1_i4, jroot) + 1_i4
    2395           0 :        a1 = -32_i4 + j1 * 16_i4
    2396             : 
    2397           0 :        j2 = (j-1_i4) / jroot
    2398           0 :        a2 = -32_i4 + j2 * 16_i4
    2399             : 
    2400           0 :        fj = real(j,dp) + (x(1) - real(a1,dp))**6 + (x(2) - real(a2,dp))**6
    2401             : 
    2402           0 :        fi = fi + 1.0_dp / fj
    2403             : 
    2404             :     end do
    2405             : 
    2406           0 :     de_jong_f5 = 1.0_dp / fi
    2407             : 
    2408           0 :     de_jong_f5 = de_jong_f5 - 0.0019996667_dp
    2409             : 
    2410           0 :   end function de_jong_f5
    2411             : 
    2412             :   ! ------------------------------------------------------------------
    2413             :   !
    2414             :   ! The Schaffer Function F6, N = 2.
    2415             :   ! Solution: x(1:n) = (/ 0.0_dp, 0.0_dp /)
    2416             :   !
    2417             :   !  Discussion:
    2418             :   !
    2419             :   !    F can be regarded as a function of R = SQRT ( X(1)^2 + X(2)^2 ).
    2420             :   !
    2421             :   !  Licensing:
    2422             :   !
    2423             :   !    This code is distributed under the GNU LGPL license.
    2424             :   !
    2425             :   !  Modified:
    2426             :   !
    2427             :   !    18 January 2001
    2428             :   !
    2429             :   !  Author:
    2430             :   !
    2431             :   !    John Burkardt
    2432             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2433             :   !
    2434             :   !  Reference:
    2435             :   !
    2436             :   !    Zbigniew Michalewicz,
    2437             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2438             :   !    Third Edition,
    2439             :   !    Springer Verlag, 1996,
    2440             :   !    ISBN: 3-540-60676-9,
    2441             :   !    LC: QA76.618.M53.
    2442             :   !
    2443             :   !  Parameters:
    2444             :   !
    2445             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2446             :   !
    2447             : 
    2448           0 :   function schaffer_f6(x)
    2449             : 
    2450             :     implicit none
    2451             : 
    2452             :     !    integer(i4) :: n
    2453             : 
    2454           0 :     real(dp) :: a
    2455           0 :     real(dp) :: b
    2456             :     real(dp) :: schaffer_f6
    2457           0 :     real(dp) :: r
    2458             :     real(dp), dimension(:), intent(in) :: x
    2459             : 
    2460           0 :     r = sqrt ( x(1)**2 + x(2)**2 )
    2461             : 
    2462           0 :     a = ( 1.0_dp + 0.001_dp * r**2 )**( -2 )
    2463             : 
    2464           0 :     b = ( sin ( r ) )**2 - 0.5_dp
    2465             : 
    2466           0 :     schaffer_f6 = 0.5_dp + a * b
    2467             : 
    2468           0 :   end function schaffer_f6
    2469             : 
    2470             :   ! ------------------------------------------------------------------
    2471             :   !
    2472             :   ! The Schaffer Function F7, N = 2.
    2473             :   ! Solution: x(1:n) = (/ 0.0_dp, 0.0_dp /)
    2474             :   !
    2475             :   !  Discussion:
    2476             :   !
    2477             :   !    Note that F is a function of R^2 = X(1)^2 + X(2)^2
    2478             :   !
    2479             :   !  Licensing:
    2480             :   !
    2481             :   !    This code is distributed under the GNU LGPL license.
    2482             :   !
    2483             :   !  Modified:
    2484             :   !
    2485             :   !    08 January 2001
    2486             :   !
    2487             :   !  Author:
    2488             :   !
    2489             :   !    John Burkardt
    2490             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2491             :   !
    2492             :   !  Reference:
    2493             :   !
    2494             :   !    Zbigniew Michalewicz,
    2495             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2496             :   !    Third Edition,
    2497             :   !    Springer Verlag, 1996,
    2498             :   !    ISBN: 3-540-60676-9,
    2499             :   !    LC: QA76.618.M53.
    2500             :   !
    2501             :   !  Parameters:
    2502             :   !
    2503             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2504             :   !
    2505             : 
    2506           0 :   function schaffer_f7(x)
    2507             : 
    2508             :     implicit none
    2509             : 
    2510             :     !    integer(i4) :: n
    2511             : 
    2512             :     real(dp) :: schaffer_f7
    2513           0 :     real(dp) :: r
    2514             :     real(dp), dimension(:), intent(in) :: x
    2515             : 
    2516           0 :     r = sqrt ( x(1)**2 + x(2)**2 )
    2517             : 
    2518           0 :     schaffer_f7 = sqrt ( r ) * ( 1.0_dp + ( sin ( 50.0_dp * r**0.2_dp ) )**2 )
    2519             : 
    2520           0 :   end function schaffer_f7
    2521             : 
    2522             :   ! ------------------------------------------------------------------
    2523             :   !
    2524             :   ! The Goldstein Price Polynomial, N = 2.
    2525             :   ! Solution:
    2526             :   !    x(1:n) = (/ 0.0_dp, -1.0_dp /)
    2527             :   !    f(x)   = 3.0_dp
    2528             :   !    http://www.pg.gda.pl/~mkwies/dyd/geadocu/fcngold.html
    2529             :   !
    2530             :   !  Discussion:
    2531             :   !
    2532             :   !    Note that F is a polynomial in X.
    2533             :   !
    2534             :   !  Licensing:
    2535             :   !
    2536             :   !    This code is distributed under the GNU LGPL license.
    2537             :   !
    2538             :   !  Modified:
    2539             :   !
    2540             :   !    08 January 2001
    2541             :   !
    2542             :   !  Author:
    2543             :   !
    2544             :   !    John Burkardt
    2545             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2546             :   !
    2547             :   !  Reference:
    2548             :   !
    2549             :   !    Zbigniew Michalewicz,
    2550             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2551             :   !    Third Edition,
    2552             :   !    Springer Verlag, 1996,
    2553             :   !    ISBN: 3-540-60676-9,
    2554             :   !    LC: QA76.618.M53.
    2555             :   !
    2556             :   !  Parameters:
    2557             :   !
    2558             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2559             :   !
    2560             : 
    2561           0 :   function goldstein_price_polynomial(x)
    2562             : 
    2563             :     implicit none
    2564             : 
    2565             :     !    integer(i4) :: n
    2566             : 
    2567           0 :     real(dp) :: a
    2568           0 :     real(dp) :: b
    2569           0 :     real(dp) :: c
    2570           0 :     real(dp) :: d
    2571             :     real(dp) :: goldstein_price_polynomial
    2572             :     real(dp), dimension(:), intent(in) :: x
    2573             : 
    2574           0 :     a = x(1) + x(2) + 1.0_dp
    2575             : 
    2576             :     b = 19.0_dp - 14.0_dp * x(1) + 3.0_dp * x(1) * x(1) - 14.0_dp * x(2) &
    2577           0 :          + 6.0_dp * x(1) * x(2) + 3.0_dp * x(2) * x(2)
    2578             : 
    2579           0 :     c = 2.0_dp * x(1) - 3.0_dp * x(2)
    2580             : 
    2581             :     d = 18.0_dp - 32.0_dp * x(1) + 12.0_dp * x(1) * x(1) + 48.0_dp * x(2) &
    2582           0 :          - 36.0_dp * x(1) * x(2) + 27.0_dp * x(2) * x(2)
    2583             : 
    2584           0 :     goldstein_price_polynomial = ( 1.0_dp + a * a * b ) * ( 30.0_dp + c * c * d )
    2585             : 
    2586           0 :     goldstein_price_polynomial = goldstein_price_polynomial - 3.0_dp
    2587             : 
    2588           0 :   end function goldstein_price_polynomial
    2589             : 
    2590             :   ! ------------------------------------------------------------------
    2591             :   !
    2592             :   ! The Branin RCOS Function, N = 2.
    2593             :   ! Solution: 1st solution: x(1:n) = (/ -pi, 12.275_dp /)
    2594             :   !           2nd solution: x(1:n) = (/  pi,  2.275_dp /)
    2595             :   !           3rd solution: x(1:n) = (/ 9.42478_dp, 2.475_dp /)
    2596             :   !
    2597             :   !  Licensing:
    2598             :   !
    2599             :   !    This code is distributed under the GNU LGPL license.
    2600             :   !
    2601             :   !  Modified:
    2602             :   !
    2603             :   !    10 January 2001
    2604             :   !
    2605             :   !  Author:
    2606             :   !
    2607             :   !    John Burkardt
    2608             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2609             :   !
    2610             :   !  Reference:
    2611             :   !
    2612             :   !    Zbigniew Michalewicz,
    2613             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2614             :   !    Third Edition,
    2615             :   !    Springer Verlag, 1996,
    2616             :   !    ISBN: 3-540-60676-9,
    2617             :   !    LC: QA76.618.M53.
    2618             :   !
    2619             :   !  Parameters:
    2620             :   !
    2621             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2622             :   !
    2623             : 
    2624           0 :   function branin_rcos(x)
    2625             : 
    2626           0 :     use mo_constants, only: pi_dp
    2627             : 
    2628             :     implicit none
    2629             : 
    2630             :     !    integer(i4) :: n
    2631             : 
    2632             :     real(dp), parameter :: a = 1.0_dp
    2633             :     real(dp) :: b
    2634             :     real(dp) :: c
    2635             :     real(dp), parameter :: d = 6.0_dp
    2636             :     real(dp), parameter :: e = 10.0_dp
    2637             :     real(dp) :: branin_rcos
    2638             :     real(dp) :: ff
    2639             :     real(dp), dimension(:), intent(in) :: x
    2640             : 
    2641           0 :     b = 5.1_dp / ( 4.0_dp * pi_dp**2 )
    2642           0 :     c = 5.0_dp / pi_dp
    2643           0 :     ff = 1.0_dp / ( 8.0_dp * pi_dp )
    2644             : 
    2645           0 :     branin_rcos = a * ( x(2) - b * x(1)**2 + c * x(1) - d )**2 &
    2646           0 :          + e * ( 1.0_dp - ff ) * cos ( x(1) ) + e
    2647             : 
    2648           0 :     branin_rcos = branin_rcos - 0.3978873577_dp
    2649             : 
    2650           0 :   end function branin_rcos
    2651             : 
    2652             :   ! ------------------------------------------------------------------
    2653             :   !
    2654             :   ! The Shekel SQRN5 Function, N = 4.
    2655             :   ! Solution: x(1:n) = (/ 4.0000371429_dp, 4.0001315700_dp, 4.0000379073_dp, 4.0001323857_dp /)
    2656             :   !
    2657             :   !  Discussion:
    2658             :   !
    2659             :   !    The minimal function value is -10.1527236935_dp.
    2660             :   !
    2661             :   !  Licensing:
    2662             :   !
    2663             :   !    This code is distributed under the GNU LGPL license.
    2664             :   !
    2665             :   !  Modified:
    2666             :   !
    2667             :   !    10 January 2001
    2668             :   !
    2669             :   !  Author:
    2670             :   !
    2671             :   !    John Burkardt
    2672             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2673             :   !
    2674             :   !  Reference:
    2675             :   !
    2676             :   !    Zbigniew Michalewicz,
    2677             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2678             :   !    Third Edition,
    2679             :   !    Springer Verlag, 1996,
    2680             :   !    ISBN: 3-540-60676-9,
    2681             :   !    LC: QA76.618.M53.
    2682             :   !
    2683             :   !  Parameters:
    2684             :   !
    2685             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2686             :   !
    2687             : 
    2688           0 :   function shekel_sqrn5(x)
    2689             : 
    2690             :     implicit none
    2691             : 
    2692             :     integer(i4), parameter :: m = 5
    2693             :     integer(i4) :: n
    2694             : 
    2695             :     real(dp), parameter, dimension ( 4, m ) :: a = reshape ( &
    2696             :          (/ 4.0_dp, 4.0_dp, 4.0_dp, 4.0_dp, &
    2697             :          1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, &
    2698             :          8.0_dp, 8.0_dp, 8.0_dp, 8.0_dp, &
    2699             :          6.0_dp, 6.0_dp, 6.0_dp, 6.0_dp, &
    2700             :          3.0_dp, 7.0_dp, 3.0_dp, 7.0_dp /), (/ 4, m /) )
    2701             :     real(dp), save, dimension ( m ) :: c = &
    2702             :          (/ 0.1_dp, 0.2_dp, 0.2_dp, 0.4_dp, 0.6_dp /)
    2703             :     real(dp) :: shekel_sqrn5
    2704             :     integer(i4) ::j
    2705             :     real(dp), dimension(:), intent(in) :: x
    2706             : 
    2707           0 :     n = size(x)
    2708           0 :     shekel_sqrn5 = 0.0_dp
    2709           0 :     do j = 1, m
    2710           0 :        shekel_sqrn5 = shekel_sqrn5 - 1.0_dp / ( c(j) + sum ( ( x(1:n) - a(1:n,j) )**2 ) )
    2711             :     end do
    2712             : 
    2713           0 :     shekel_sqrn5 = shekel_sqrn5 + 10.1527236935_dp
    2714             : 
    2715           0 :   end function shekel_sqrn5
    2716             : 
    2717             :   ! ------------------------------------------------------------------
    2718             :   !
    2719             :   ! The Shekel SQRN7 Function, N = 4.
    2720             :   ! Solution: x(1:n) = (/ 4.0005729560_dp, 4.0006881764_dp, 3.9994902225_dp, 3.9996048794_dp /)
    2721             :   !
    2722             :   !  Licensing:
    2723             :   !
    2724             :   !    This code is distributed under the GNU LGPL license.
    2725             :   !
    2726             :   !  Modified:
    2727             :   !
    2728             :   !    12 January 2001
    2729             :   !
    2730             :   !  Author:
    2731             :   !
    2732             :   !    John Burkardt
    2733             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2734             :   !
    2735             :   !  Reference:
    2736             :   !
    2737             :   !    Zbigniew Michalewicz,
    2738             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2739             :   !    Third Edition,
    2740             :   !    Springer Verlag, 1996,
    2741             :   !    ISBN: 3-540-60676-9,
    2742             :   !    LC: QA76.618.M53.
    2743             :   !
    2744             :   !  Parameters:
    2745             :   !
    2746             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2747             :   !
    2748             : 
    2749           0 :   function shekel_sqrn7(x)
    2750             : 
    2751             :     implicit none
    2752             : 
    2753             :     integer(i4), parameter :: m = 7
    2754             :     integer(i4) :: n
    2755             : 
    2756             :     real(dp), parameter, dimension ( 4, m ) :: a = reshape ( &
    2757             :          (/ 4.0_dp, 4.0_dp, 4.0_dp, 4.0_dp, &
    2758             :          1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, &
    2759             :          8.0_dp, 8.0_dp, 8.0_dp, 8.0_dp, &
    2760             :          6.0_dp, 6.0_dp, 6.0_dp, 6.0_dp, &
    2761             :          3.0_dp, 7.0_dp, 3.0_dp, 7.0_dp, &
    2762             :          2.0_dp, 9.0_dp, 2.0_dp, 9.0_dp, &
    2763             :          5.0_dp, 5.0_dp, 3.0_dp, 3.0_dp /), (/ 4, m /) )
    2764             :     real(dp), save, dimension ( m ) :: c = &
    2765             :          (/ 0.1_dp, 0.2_dp, 0.2_dp, 0.4_dp, 0.6_dp, 0.6_dp, 0.3_dp /)
    2766             :     real(dp) :: shekel_sqrn7
    2767             :     integer(i4) ::j
    2768             :     real(dp), dimension(:), intent(in) :: x
    2769             : 
    2770           0 :     n = size(x)
    2771           0 :     shekel_sqrn7 = 0.0_dp
    2772           0 :     do j = 1, m
    2773           0 :        shekel_sqrn7 = shekel_sqrn7 - 1.0_dp / ( c(j) + sum ( ( x(1:n) - a(1:n,j) )**2 ) )
    2774             :     end do
    2775             : 
    2776           0 :     shekel_sqrn7 = shekel_sqrn7 + 10.4024645722_dp
    2777             : 
    2778           0 :   end function shekel_sqrn7
    2779             : 
    2780             :   ! ------------------------------------------------------------------
    2781             :   !
    2782             :   ! The Shekel SQRN10 Function, N = 4.
    2783             :   ! Solution: x(1:n) = (/ 4.0007465727_dp, 4.0005916919_dp, 3.9996634360_dp, 3.9995095935_dp /)
    2784             :   !
    2785             :   !  Licensing:
    2786             :   !
    2787             :   !    This code is distributed under the GNU LGPL license.
    2788             :   !
    2789             :   !  Modified:
    2790             :   !
    2791             :   !    12 January 2001
    2792             :   !
    2793             :   !  Author:
    2794             :   !
    2795             :   !    John Burkardt
    2796             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2797             :   !
    2798             :   !  Reference:
    2799             :   !
    2800             :   !    Zbigniew Michalewicz,
    2801             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2802             :   !    Third Edition,
    2803             :   !    Springer Verlag, 1996,
    2804             :   !    ISBN: 3-540-60676-9,
    2805             :   !    LC: QA76.618.M53.
    2806             :   !
    2807             :   !  Parameters:
    2808             :   !
    2809             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2810             :   !
    2811             : 
    2812           0 :   function shekel_sqrn10(x)
    2813             : 
    2814             :     implicit none
    2815             : 
    2816             :     integer(i4), parameter :: m = 10
    2817             :     integer(i4) :: n
    2818             : 
    2819             :     real(dp), parameter, dimension ( 4, m ) :: a = reshape ( &
    2820             :          (/ 4.0, 4.0, 4.0, 4.0, &
    2821             :          1.0, 1.0, 1.0, 1.0, &
    2822             :          8.0, 8.0, 8.0, 8.0, &
    2823             :          6.0, 6.0, 6.0, 6.0, &
    2824             :          3.0, 7.0, 3.0, 7.0, &
    2825             :          2.0, 9.0, 2.0, 9.0, &
    2826             :          5.0, 5.0, 3.0, 3.0, &
    2827             :          8.0, 1.0, 8.0, 1.0, &
    2828             :          6.0, 2.0, 6.0, 2.0, &
    2829             :          7.0, 3.6, 7.0, 3.6 /), (/ 4, m /) )
    2830             : 
    2831             :     real(dp), save, dimension ( m ) :: c = &
    2832             :          (/ 0.1_dp, 0.2_dp, 0.2_dp, 0.4_dp, 0.6_dp, &
    2833             :          0.6_dp, 0.3_dp, 0.7_dp, 0.5_dp, 0.5_dp /)
    2834             :     real(dp) :: shekel_sqrn10
    2835             :     integer(i4) ::j
    2836             :     real(dp), dimension(:), intent(in) :: x
    2837             : 
    2838           0 :     n = size(x)
    2839           0 :     shekel_sqrn10 = 0.0_dp
    2840           0 :     do j = 1, m
    2841           0 :        shekel_sqrn10 = shekel_sqrn10 - 1.0_dp / ( c(j) + sum ( ( x(1:n) - a(1:n,j) )**2 ) )
    2842             :     end do
    2843             : 
    2844           0 :     shekel_sqrn10 = shekel_sqrn10 + 10.5359339075_dp
    2845             : 
    2846           0 :   end function shekel_sqrn10
    2847             : 
    2848             :   ! ------------------------------------------------------------------
    2849             :   !
    2850             :   ! The Six-Hump Camel-Back Polynomial, N = 2.
    2851             :   ! Solution: 1st solution: x(1:n) = (/ -0.0898_dp,  0.7126_dp /)
    2852             :   !           2nd solution: x(1:n) = (/  0.0898_dp, -0.7126_dp /)
    2853             :   !
    2854             :   !  Licensing:
    2855             :   !
    2856             :   !    This code is distributed under the GNU LGPL license.
    2857             :   !
    2858             :   !  Modified:
    2859             :   !
    2860             :   !    12 January 2001
    2861             :   !
    2862             :   !  Author:
    2863             :   !
    2864             :   !    John Burkardt
    2865             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2866             :   !
    2867             :   !  Reference:
    2868             :   !
    2869             :   !    Zbigniew Michalewicz,
    2870             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2871             :   !    Third Edition,
    2872             :   !    Springer Verlag, 1996,
    2873             :   !    ISBN: 3-540-60676-9,
    2874             :   !    LC: QA76.618.M53.
    2875             :   !
    2876             :   !  Parameters:
    2877             :   !
    2878             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2879             :   !
    2880             : 
    2881           0 :   function six_hump_camel_back_polynomial(x)
    2882             : 
    2883             :     implicit none
    2884             : 
    2885             :     !    integer(i4) :: n
    2886             : 
    2887             :     real(dp) :: six_hump_camel_back_polynomial
    2888             :     real(dp), dimension(:), intent(in) :: x
    2889             : 
    2890           0 :     six_hump_camel_back_polynomial = ( 4.0_dp - 2.1_dp * x(1)**2 + x(1)**4 / 3.0_dp ) * x(1)**2 &
    2891           0 :          + x(1) * x(2) + 4.0_dp * ( x(2)**2 - 1.0_dp ) * x(2)**2
    2892             : 
    2893           0 :     six_hump_camel_back_polynomial = six_hump_camel_back_polynomial + 1.0316284229_dp
    2894             : 
    2895           0 :   end function six_hump_camel_back_polynomial
    2896             : 
    2897             :   ! ------------------------------------------------------------------
    2898             :   !
    2899             :   ! The Shubert Function, N = 2.
    2900             :   ! Solution: x(1:n) = (/ -7.0835064076515595_dp, -7.708313735499347_dp /)
    2901             :   !
    2902             :   !  Discussion:
    2903             :   !
    2904             :   !    For -10 <= X(I) <= 10, the function has 760 local minima, 18 of which
    2905             :   !    are global minima, with minimum value -186.73090883102378.
    2906             :   !
    2907             :   !  Licensing:
    2908             :   !
    2909             :   !    This code is distributed under the GNU LGPL license.
    2910             :   !
    2911             :   !  Modified:
    2912             :   !
    2913             :   !    12 January 2001
    2914             :   !
    2915             :   !  Author:
    2916             :   !
    2917             :   !    John Burkardt
    2918             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2919             :   !
    2920             :   !  Reference:
    2921             :   !
    2922             :   !    Zbigniew Michalewicz,
    2923             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2924             :   !    Third Edition,
    2925             :   !    Springer Verlag, 1996,
    2926             :   !    ISBN: 3-540-60676-9,
    2927             :   !    LC: QA76.618.M53.
    2928             :   !
    2929             :   !    Bruno Shubert,
    2930             :   !    A sequential method seeking the global maximum of a function,
    2931             :   !    SIAM Journal on Numerical Analysis,
    2932             :   !    Volume 9, pages 379-388, 1972.
    2933             :   !
    2934             :   !  Parameters:
    2935             :   !
    2936             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2937             :   !
    2938             : 
    2939           0 :   function shubert(x)
    2940             : 
    2941             :     implicit none
    2942             : 
    2943             :     integer(i4) :: n
    2944             : 
    2945             :     real(dp) :: shubert
    2946           0 :     real(dp) :: factor
    2947             :     integer(i4) ::i
    2948             :     integer(i4) ::k
    2949           0 :     real(dp) :: k_r8
    2950             :     real(dp), dimension(:), intent(in) :: x
    2951             : 
    2952           0 :     n = size(x)
    2953           0 :     shubert = 1.0_dp
    2954             : 
    2955           0 :     do i = 1, n
    2956             :        factor = 0.0_dp
    2957           0 :        do k = 1, 5
    2958           0 :           k_r8 = real ( k, dp )
    2959           0 :           factor = factor + k_r8 * cos ( ( k_r8 + 1.0_dp ) * x(i) + k_r8 )
    2960             :        end do
    2961           0 :        shubert = shubert * factor
    2962             :     end do
    2963             : 
    2964           0 :     shubert = shubert + 186.7309088310_dp
    2965             : 
    2966           0 :   end function shubert
    2967             : 
    2968             :   ! ------------------------------------------------------------------
    2969             :   !
    2970             :   ! The Stuckman Function, N = 2.
    2971             :   ! Solution: Only iterative solution; Check reference.
    2972             :   !
    2973             :   !  Licensing:
    2974             :   !
    2975             :   !    This code is distributed under the GNU LGPL license.
    2976             :   !
    2977             :   !  Modified:
    2978             :   !
    2979             :   !    16 October 2004
    2980             :   !
    2981             :   !  Author:
    2982             :   !
    2983             :   !    John Burkardt
    2984             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    2985             :   !
    2986             :   !  Reference:
    2987             :   !
    2988             :   !    Zbigniew Michalewicz,
    2989             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    2990             :   !    Third Edition,
    2991             :   !    Springer Verlag, 1996,
    2992             :   !    ISBN: 3-540-60676-9,
    2993             :   !    LC: QA76.618.M53.
    2994             :   !
    2995             :   !  Parameters:
    2996             :   !
    2997             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    2998             :   !
    2999             : 
    3000           0 :   function stuckman(x)
    3001             : 
    3002           0 :     use mo_utils, only: eq, le
    3003             : 
    3004             :     implicit none
    3005             : 
    3006             :     !    integer(i4) :: n
    3007             : 
    3008           0 :     real(dp) :: a1
    3009           0 :     real(dp) :: a2
    3010           0 :     real(dp) :: b
    3011             :     real(dp) :: stuckman
    3012           0 :     real(dp) :: m1
    3013           0 :     real(dp) :: m2
    3014           0 :     real(dp) :: r11
    3015           0 :     real(dp) :: r12
    3016           0 :     real(dp) :: r21
    3017           0 :     real(dp) :: r22
    3018             :     real(dp), dimension(:), intent(in) :: x
    3019             : 
    3020             :     real(dp), save :: b_save = 0.0_dp
    3021             :     real(dp), save :: m1_save = 0.0_dp
    3022             :     real(dp), save :: m2_save = 0.0_dp
    3023             :     real(dp), save :: r11_save = 0.0_dp
    3024             :     real(dp), save :: r12_save = 0.0_dp
    3025             :     real(dp), save :: r21_save = 0.0_dp
    3026             :     real(dp), save :: r22_save = 0.0_dp
    3027             : 
    3028             :     !call p36_p_get ( b, m1, m2, r11, r12, r21, r22, seed )
    3029           0 :     b = b_save
    3030           0 :     m1 = m1_save
    3031           0 :     m2 = m2_save
    3032           0 :     r11 = r11_save
    3033           0 :     r12 = r12_save
    3034           0 :     r21 = r21_save
    3035           0 :     r22 = r22_save
    3036             : 
    3037           0 :     a1 = r8_aint ( abs ( x(1) - r11 ) ) + r8_aint ( abs ( x(2) - r21 ) )
    3038           0 :     a2 = r8_aint ( abs ( x(1) - r12 ) ) + r8_aint ( abs ( x(2) - r22 ) )
    3039             : 
    3040           0 :     if ( le(x(1),b) ) then
    3041           0 :        if ( eq(a1,0.0_dp) ) then
    3042           0 :           stuckman = r8_aint ( m1 )
    3043             :        else
    3044           0 :           stuckman = r8_aint ( m1 * sin ( a1 ) / a1 )
    3045             :        end if
    3046             :     else
    3047           0 :        if ( eq(a2,0.0_dp) ) then
    3048           0 :           stuckman = r8_aint ( m2 )
    3049             :        else
    3050           0 :           stuckman = r8_aint ( m2 * sin ( a2 ) / a2 )
    3051             :        end if
    3052             :     end if
    3053             : 
    3054           0 :   end function stuckman
    3055             : 
    3056             :   ! ------------------------------------------------------------------
    3057             :   !
    3058             :   ! The Easom Function, N = 2.
    3059             :   ! Solution: x(1:n) = (/ pi, pi /)
    3060             :   !
    3061             :   !  Licensing:
    3062             :   !
    3063             :   !    This code is distributed under the GNU LGPL license.
    3064             :   !
    3065             :   !  Modified:
    3066             :   !
    3067             :   !    11 January 2001
    3068             :   !
    3069             :   !  Author:
    3070             :   !
    3071             :   !    John Burkardt
    3072             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3073             :   !
    3074             :   !  Reference:
    3075             :   !
    3076             :   !    Zbigniew Michalewicz,
    3077             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    3078             :   !    Third Edition,
    3079             :   !    Springer Verlag, 1996,
    3080             :   !    ISBN: 3-540-60676-9,
    3081             :   !    LC: QA76.618.M53.
    3082             :   !
    3083             :   !  Parameters:
    3084             :   !
    3085             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3086             :   !
    3087             : 
    3088           0 :   function easom(x)
    3089             : 
    3090           0 :     use mo_constants, only: pi_dp
    3091             : 
    3092             :     implicit none
    3093             : 
    3094             :     !    integer(i4) :: n
    3095             : 
    3096           0 :     real(dp) :: arg
    3097             :     real(dp) :: easom
    3098             :     real(dp), dimension(:), intent(in) :: x
    3099             : 
    3100           0 :     arg = - ( x(1) - pi_dp )**2 - ( x(2) - pi_dp )**2
    3101           0 :     easom = - cos ( x(1) ) * cos ( x(2) ) * exp ( arg )
    3102           0 :     easom = easom + 1.0_dp
    3103             : 
    3104           0 :   end function easom
    3105             : 
    3106             :   ! ------------------------------------------------------------------
    3107             :   !
    3108             :   ! The Bohachevsky Function #1, N = 2.
    3109             :   ! Solution: x(1:n) = (/ 0.0_dp, 0.0_dp /)
    3110             :   !
    3111             :   !  Licensing:
    3112             :   !
    3113             :   !    This code is distributed under the GNU LGPL license.
    3114             :   !
    3115             :   !  Modified:
    3116             :   !
    3117             :   !    11 January 2001
    3118             :   !
    3119             :   !  Author:
    3120             :   !
    3121             :   !    John Burkardt
    3122             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3123             :   !
    3124             :   !  Reference:
    3125             :   !
    3126             :   !    Zbigniew Michalewicz,
    3127             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    3128             :   !    Third Edition,
    3129             :   !    Springer Verlag, 1996,
    3130             :   !    ISBN: 3-540-60676-9,
    3131             :   !    LC: QA76.618.M53.
    3132             :   !
    3133             :   !  Parameters:
    3134             :   !
    3135             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3136             :   !
    3137             : 
    3138           0 :   function bohachevsky1(x)
    3139             : 
    3140           0 :     use mo_constants, only: pi_dp
    3141             : 
    3142             :     implicit none
    3143             : 
    3144             :     !    integer(i4) :: n
    3145             : 
    3146             :     real(dp) :: bohachevsky1
    3147             :     real(dp), dimension(:), intent(in) :: x
    3148             : 
    3149           0 :     bohachevsky1 =           x(1) * x(1) - 0.3_dp * cos ( 3.0_dp * pi_dp * x(1) ) &
    3150           0 :          + 2.0_dp * x(2) * x(2) - 0.4_dp * cos ( 4.0_dp * pi_dp * x(2) ) &
    3151           0 :          + 0.7_dp
    3152             : 
    3153           0 :   end function bohachevsky1
    3154             : 
    3155             :   ! ------------------------------------------------------------------
    3156             :   !
    3157             :   ! The Bohachevsky Function #2, N = 2.
    3158             :   ! Solution: x(1:n) = (/ 0.0_dp, 0.0_dp /)
    3159             :   !
    3160             :   !  Licensing:
    3161             :   !
    3162             :   !    This code is distributed under the GNU LGPL license.
    3163             :   !
    3164             :   !  Modified:
    3165             :   !
    3166             :   !    11 January 2001
    3167             :   !
    3168             :   !  Author:
    3169             :   !
    3170             :   !    John Burkardt
    3171             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3172             :   !
    3173             :   !  Reference:
    3174             :   !
    3175             :   !    Zbigniew Michalewicz,
    3176             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    3177             :   !    Third Edition,
    3178             :   !    Springer Verlag, 1996,
    3179             :   !    ISBN: 3-540-60676-9,
    3180             :   !    LC: QA76.618.M53.
    3181             :   !
    3182             :   !  Parameters:
    3183             :   !
    3184             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3185             :   !
    3186             : 
    3187           0 :   function bohachevsky2(x)
    3188             : 
    3189           0 :     use mo_constants, only: pi_dp
    3190             : 
    3191             :     implicit none
    3192             : 
    3193             :     !    integer(i4) :: n
    3194             : 
    3195             :     real(dp) :: bohachevsky2
    3196             :     real(dp), dimension(:), intent(in) :: x
    3197             : 
    3198           0 :     bohachevsky2 = x(1) * x(1) + 2.0_dp * x(2) * x(2) &
    3199             :          - 0.3_dp * cos ( 3.0_dp * pi_dp * x(1) ) &
    3200           0 :          * cos ( 4.0_dp * pi_dp * x(2) ) + 0.3_dp
    3201             : 
    3202           0 :   end function bohachevsky2
    3203             : 
    3204             :   ! ------------------------------------------------------------------
    3205             :   !
    3206             :   ! The Bohachevsky Function #3, N = 2.
    3207             :   ! Solution:
    3208             :   !     x(1:2) = (/ 0.0_dp, 0.0_dp /)
    3209             :   !     f(x)   = 0.0_dp
    3210             :   !
    3211             :   !  Discussion:
    3212             :   !    J. Burkhardt:
    3213             :   !       There is a typo in the reference.  I'm just guessing at the correction.
    3214             :   !    J. Mai:
    3215             :   !       Typo in function found.
    3216             :   !       see: http://www-optima.amp.i.kyoto-u.ac.jp/member/student/
    3217             :   !            hedar/Hedar_files/TestGO_files/Page595.htm
    3218             :   !
    3219             :   !  Licensing:
    3220             :   !
    3221             :   !    This code is distributed under the GNU LGPL license.
    3222             :   !
    3223             :   !  Modified:
    3224             :   !
    3225             :   !    12 January 2001
    3226             :   !
    3227             :   !  Author:
    3228             :   !
    3229             :   !    John Burkardt
    3230             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3231             :   !
    3232             :   !  Reference:
    3233             :   !
    3234             :   !    Zbigniew Michalewicz,
    3235             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    3236             :   !    Third Edition,
    3237             :   !    Springer Verlag, 1996,
    3238             :   !    ISBN: 3-540-60676-9,
    3239             :   !    LC: QA76.618.M53.
    3240             :   !
    3241             :   !  Parameters:
    3242             :   !
    3243             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3244             :   !
    3245             : 
    3246           0 :   function bohachevsky3(x)
    3247             : 
    3248           0 :     use mo_constants, only: pi_dp
    3249             : 
    3250             :     implicit none
    3251             : 
    3252             :     !    integer(i4) :: n
    3253             : 
    3254             :     real(dp) :: bohachevsky3
    3255             :     real(dp), dimension(:), intent(in) :: x
    3256             : 
    3257           0 :     bohachevsky3 = x(1)**2 + 2.0_dp * x(2)**2 &
    3258             :          - 0.3_dp * cos ( 3.0_dp * pi_dp * x(1)  &
    3259           0 :          + 4.0_dp * pi_dp * x(2) ) + 0.3_dp
    3260             : 
    3261           0 :   end function bohachevsky3
    3262             : 
    3263             :   ! ------------------------------------------------------------------
    3264             :   !
    3265             :   ! The Colville Polynomial, N = 4.
    3266             :   ! Solution: x(1:n) = (/ 1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp /)
    3267             :   !
    3268             :   !  Licensing:
    3269             :   !
    3270             :   !    This code is distributed under the GNU LGPL license.
    3271             :   !
    3272             :   !  Modified:
    3273             :   !
    3274             :   !    12 January 2001
    3275             :   !
    3276             :   !  Author:
    3277             :   !
    3278             :   !    John Burkardt
    3279             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3280             :   !
    3281             :   !  Reference:
    3282             :   !
    3283             :   !    Zbigniew Michalewicz,
    3284             :   !    Genetic Algorithms + Data Structures = Evolution Programs,
    3285             :   !    Third Edition,
    3286             :   !    Springer Verlag, 1996,
    3287             :   !    ISBN: 3-540-60676-9,
    3288             :   !    LC: QA76.618.M53.
    3289             :   !
    3290             :   !  Parameters:
    3291             :   !
    3292             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3293             :   !
    3294             : 
    3295           0 :   function colville_polynomial(x)
    3296             : 
    3297             :     implicit none
    3298             : 
    3299             :     !    integer(i4) :: n
    3300             : 
    3301             :     real(dp) :: colville_polynomial
    3302             :     real(dp), dimension(:), intent(in) :: x
    3303             : 
    3304           0 :     colville_polynomial = 100.0_dp * ( x(2) - x(1)**2 )**2 &
    3305           0 :          + ( 1.0_dp - x(1) )**2 &
    3306           0 :          + 90.0_dp * ( x(4) - x(3)**2 )**2 &
    3307           0 :          + ( 1.0_dp - x(3) )**2 &
    3308             :          + 10.1_dp * ( ( x(2) - 1.0_dp )**2 + ( x(4) - 1.0_dp )**2 ) &
    3309           0 :          + 19.8_dp * ( x(2) - 1.0_dp ) * ( x(4) - 1.0_dp )
    3310             : 
    3311           0 :   end function colville_polynomial
    3312             : 
    3313             :   ! ------------------------------------------------------------------
    3314             :   !
    3315             :   ! The Powell 3D function, N = 3.
    3316             :   ! Solution: x(1:n) = (/ 1.0_dp, 1.0_dp, 1.0_dp /)
    3317             :   !
    3318             :   !  Licensing:
    3319             :   !
    3320             :   !    This code is distributed under the GNU LGPL license.
    3321             :   !
    3322             :   !  Modified:
    3323             :   !
    3324             :   !    03 March 2002
    3325             :   !
    3326             :   !  Author:
    3327             :   !
    3328             :   !    John Burkardt
    3329             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3330             :   !
    3331             :   !  Reference:
    3332             :   !
    3333             :   !    MJD Powell,
    3334             :   !    An Efficient Method for Finding the Minimum of a Function of
    3335             :   !    Several Variables Without Calculating Derivatives,
    3336             :   !    Computer Journal,
    3337             :   !    Volume 7, Number 2, pages 155-162, 1964.
    3338             :   !
    3339             :   !  Parameters:
    3340             :   !
    3341             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3342             :   !
    3343             : 
    3344           0 :   function powell3d(x)
    3345             : 
    3346           0 :     use mo_constants, only: pi_dp
    3347             :     use mo_utils, only: eq
    3348             : 
    3349             :     implicit none
    3350             : 
    3351             :     !    integer(i4) :: n
    3352             : 
    3353           0 :     real(dp) :: arg
    3354             :     real(dp) :: powell3d
    3355           0 :     real(dp) :: term
    3356             :     real(dp), dimension(:), intent(in) :: x
    3357             : 
    3358           0 :     if ( eq(x(2),0.0_dp) ) then
    3359             :        term = 0.0_dp
    3360             :     else
    3361             :        !arg = ( x(1) + 2.0_dp * x(2) + x(3) ) / x(2)
    3362             :        ! changed according to original paper of Powell (1964)
    3363           0 :        arg = ( x(1) + x(3) ) / x(2) - 2.0_dp
    3364             : 
    3365           0 :        term = arg*arg
    3366           0 :        if ( term .lt. 708._dp ) then
    3367           0 :           term = exp ( - term )
    3368             :        else
    3369             :           ! avoid underflow
    3370             :           term = 0.0_dp
    3371             :        end if
    3372             :     end if
    3373             : 
    3374             :     powell3d = 3.0_dp &
    3375           0 :          - 1.0_dp / ( 1.0_dp + ( x(1) - x(2) )**2 ) &
    3376           0 :          - sin ( 0.5_dp * pi_dp * x(2) * x(3) ) &
    3377           0 :          - term
    3378             : 
    3379           0 :   end function powell3d
    3380             : 
    3381             :   ! ------------------------------------------------------------------
    3382             :   !
    3383             :   ! The Himmelblau function, N = 2.
    3384             :   ! Solution: x(1:2) = (/ 3.0_dp, 2.0_dp /)
    3385             :   !
    3386             :   !  Discussion:
    3387             :   !
    3388             :   !    This function has 4 global minima:
    3389             :   !
    3390             :   !      X* = (  3,        2       ), F(X*) = 0.
    3391             :   !      X* = (  3.58439, -1.84813 ), F(X*) = 0.
    3392             :   !      X* = ( -3.77934, -3.28317 ), F(X*) = 0.
    3393             :   !      X* = ( -2.80512,  3.13134 ), F(X*) = 0.
    3394             :   !
    3395             :   !  Licensing:
    3396             :   !
    3397             :   !    This code is distributed under the GNU LGPL license.
    3398             :   !
    3399             :   !  Modified:
    3400             :   !
    3401             :   !    28 January 2008
    3402             :   !
    3403             :   !  Author:
    3404             :   !
    3405             :   !    John Burkardt
    3406             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3407             :   !
    3408             :   !  Reference:
    3409             :   !
    3410             :   !    David Himmelblau,
    3411             :   !    Applied Nonlinear Programming,
    3412             :   !    McGraw Hill, 1972,
    3413             :   !    ISBN13: 978-0070289215,
    3414             :   !   LC: T57.8.H55.
    3415             :   !
    3416             :   !  Parameters:
    3417             :   !
    3418             :   !    Input, real(dp) :: X(N), the argument of the objective function.
    3419             :   !
    3420             : 
    3421           0 :   function himmelblau(x)
    3422             : 
    3423             :     implicit none
    3424             : 
    3425             :     !    integer(i4) :: n
    3426             : 
    3427             :     real(dp) :: himmelblau
    3428             :     real(dp), dimension(:), intent(in) :: x
    3429             : 
    3430           0 :     himmelblau = ( x(1)**2 + x(2) - 11.0_dp )**2 &
    3431           0 :          + ( x(1) + x(2)**2 - 7.0_dp )**2
    3432             : 
    3433           0 :   end function himmelblau
    3434             : 
    3435             :   ! ------------------------------------------------------------------
    3436             :   !
    3437             :   ! The Griewank Function, N=2 or N=10.
    3438             :   ! Solution: x(1:n) = 0.0_dp
    3439             : 
    3440             :   !
    3441             :   ! Coded originally by Q Duan. Edited for incorporation into Fortran DDS algorithm by
    3442             :   ! Bryan Tolson, Nov 2005.
    3443             :   !
    3444             :   ! Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3445             :   !
    3446             :   ! I/O Variable definitions:
    3447             :   !     nopt     -  the number of decision variables
    3448             :   !     x_values -      an array of decision variable values (size nopt)
    3449             :   !     fvalue   -      the value of the objective function with x_values as input
    3450             : 
    3451           0 :   function griewank(x_values)
    3452             : 
    3453           0 :     use mo_kind, only: i4, dp
    3454             : 
    3455             :     implicit none
    3456             : 
    3457             :     real(dp), dimension(:), intent(in)  :: x_values
    3458             :     real(dp) :: griewank
    3459             : 
    3460             :     integer(i4) :: nopt
    3461             :     integer(i4) :: j
    3462           0 :     real(dp)    :: d, u1, u2
    3463             : 
    3464           0 :     nopt = size(x_values)
    3465           0 :     if (nopt .eq. 2) then
    3466             :        d = 200.0_dp
    3467             :     else
    3468           0 :        d = 4000.0_dp
    3469             :     end if
    3470           0 :     u1 = sum(x_values**2) / d
    3471           0 :     u2 = 1.0_dp
    3472           0 :     do j=1, nopt
    3473           0 :        u2 = u2 * cos(x_values(j)/sqrt(real(j,dp)))
    3474             :     end do
    3475           0 :     griewank = u1 - u2 + 1.0_dp
    3476             :     !
    3477           0 :   end function griewank
    3478             : 
    3479             :   ! ------------------------------------------------------------------
    3480             :   !
    3481             :   !  The Rosenbrock parabolic value function, N = 2.
    3482             :   !  Solution: x(1:n) = 1.0_dp
    3483             :   !
    3484             :   !  Licensing:
    3485             :   !
    3486             :   !    This code is distributed under the GNU LGPL license.
    3487             :   !
    3488             :   !  Modified:
    3489             :   !
    3490             :   !    19 February 2008
    3491             :   !
    3492             :   !  Author:
    3493             :   !
    3494             :   !    John Burkardt
    3495             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3496             :   !
    3497             :   !  Reference:
    3498             :   !
    3499             :   !    R ONeill,
    3500             :   !    Algorithm AS 47:
    3501             :   !    Function Minimization Using a Simplex Procedure,
    3502             :   !    Applied Statistics,
    3503             :   !    Volume 20, Number 3, 1971, pages 338-345.
    3504             :   !
    3505             :   !  Parameters:
    3506             :   !
    3507             :   !    Input, real(dp) X(2), the argument.
    3508             :   !
    3509             : 
    3510           0 :   function rosenbrock(x)
    3511             : 
    3512             :     implicit none
    3513             : 
    3514             :     real(dp), dimension(:), intent(in) :: x
    3515             :     real(dp) :: rosenbrock
    3516             : 
    3517           0 :     real(dp) :: fx
    3518           0 :     real(dp) :: fx1
    3519           0 :     real(dp) :: fx2
    3520             : 
    3521           0 :     fx1 = x(2) - x(1) * x(1)
    3522           0 :     fx2 = 1.0_dp - x(1)
    3523             : 
    3524           0 :     fx = 100.0_dp * fx1 * fx1 + fx2 * fx2
    3525             : 
    3526           0 :     rosenbrock = fx
    3527             : 
    3528           0 :   end function rosenbrock
    3529             : 
    3530             :   ! ------------------------------------------------------------------
    3531             :   !
    3532             :   !  The Sphere model, N >= 1.
    3533             :   !  Solution: x(1:n) = 0.0_dp
    3534             :   !
    3535             :   !  Author:
    3536             :   !
    3537             :   !    Matlab code by A. Hedar
    3538             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3539             :   !
    3540             :   !  Parameters:
    3541             :   !
    3542             :   !    Input, real(dp) X(N), the argument.
    3543             :   !
    3544             : 
    3545           0 :   function sphere_model(x)
    3546             : 
    3547             :     implicit none
    3548             : 
    3549             :     real(dp), dimension(:), intent(in) :: x
    3550             :     real(dp) :: sphere_model
    3551             : 
    3552             :     integer(i4) :: n
    3553             :     integer(i4) :: j
    3554             : 
    3555           0 :     n = size(x)
    3556           0 :     sphere_model = 0.0_dp
    3557           0 :     do j=1, n
    3558           0 :        sphere_model = sphere_model + x(j)**2
    3559             :     enddo
    3560             : 
    3561           0 :   end function sphere_model
    3562             : 
    3563             :   ! ------------------------------------------------------------------
    3564             :   !
    3565             :   !  The Rastrigin function, N >= 2.
    3566             :   !  Solution: x(1:n) = 0.0_dp
    3567             :   !  Search domain: -5.12 <= xi <= 5.12
    3568             :   !
    3569             :   !  Author:
    3570             :   !
    3571             :   !    Matlab code by A. Hedar
    3572             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3573             :   !
    3574             :   !  Parameters:
    3575             :   !
    3576             :   !    Input, real(dp) X(N), the argument.
    3577             :   !
    3578             : 
    3579           0 :   function rastrigin(x)
    3580             : 
    3581           0 :     use mo_constants, only: pi_dp
    3582             : 
    3583             :     implicit none
    3584             : 
    3585             :     real(dp), dimension(:), intent(in) :: x
    3586             :     real(dp) :: rastrigin
    3587             : 
    3588             :     integer(i4) :: n
    3589             :     integer(i4) :: j
    3590             : 
    3591           0 :     n = size(x,dim=1)
    3592           0 :     rastrigin = 0.0_dp
    3593           0 :     do j=1, n
    3594           0 :        rastrigin = rastrigin+ (x(j)**2 - 10.0_dp*cos(2.0_dp*pi_dp*x(j)))
    3595             :     enddo
    3596           0 :     rastrigin = rastrigin + 10.0_dp * real(n,dp)
    3597             : 
    3598             :     ! FUNCTN04 of SCEUA F77 source code
    3599             :     !   Bound: X1=[-1,1], X2=[-1,1]
    3600             :     !   Global Optimum: -2, (0,0)
    3601             :     ! n = size(x,dim=1)
    3602             :     ! rastrigin = 0.0_dp
    3603             :     ! do j=1, n
    3604             :     !    rastrigin = rastrigin+ (x(j)**2 - cos(18.0_dp*x(j)))
    3605             :     ! enddo
    3606             : 
    3607           0 :   end function rastrigin
    3608             : 
    3609             :   ! ------------------------------------------------------------------
    3610             :   !
    3611             :   !  The Schwefel function, N >= 2.
    3612             :   !  Solution: x(1:n) = 1.0_dp
    3613             :   !  Solution: x(1:n) = 420.968746_dp   ( see (2) and (3) )
    3614             :   !
    3615             :   !  Author:
    3616             :   !
    3617             :   !    (1) Matlab code by A. Hedar
    3618             :   !    (2) http://www.aridolan.com/ga/gaa/Schwefel.html
    3619             :   !    (3) http://www.it.lut.fi/ip/evo/functions/node10.html
    3620             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3621             :   !
    3622             :   !  Parameters:
    3623             :   !
    3624             :   !    Input, real(dp) X(N), the argument.
    3625             :   !
    3626             : 
    3627           0 :   function schwefel(x)
    3628             : 
    3629             :     implicit none
    3630             : 
    3631             :     real(dp), dimension(:), intent(in) :: x
    3632             :     real(dp) :: schwefel
    3633             : 
    3634             :     integer(i4) :: n
    3635             : 
    3636           0 :     n = size(x)
    3637           0 :     schwefel = sum(-x*sin(sqrt(abs(x)))) + 418.982887272433799807913601398_dp*real(n,dp)
    3638             : 
    3639           0 :   end function schwefel
    3640             : 
    3641             :   ! ------------------------------------------------------------------
    3642             :   !
    3643             :   !  The Ackley function, N >= 2.
    3644             :   !  Solution: x(1:n) = 0.0_dp
    3645             :   !
    3646             :   !  Author:
    3647             :   !
    3648             :   !    Matlab code by A. Hedar
    3649             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3650             :   !
    3651             :   !  Parameters:
    3652             :   !
    3653             :   !    Input, real(dp) X(N), the argument.
    3654             :   !
    3655             : 
    3656           0 :   function ackley(x)
    3657             : 
    3658           0 :     use mo_constants, only: pi_dp
    3659             : 
    3660             :     implicit none
    3661             : 
    3662             :     real(dp), dimension(:), intent(in) :: x
    3663             :     real(dp) :: ackley
    3664             : 
    3665             :     integer(i4) :: n
    3666             :     real(dp), parameter :: a = 20.0_dp
    3667             :     real(dp), parameter :: b = 0.2_dp
    3668             :     real(dp), parameter :: c = 2.0_dp*pi_dp
    3669           0 :     real(dp) :: s1, s2
    3670             : 
    3671           0 :     n = size(x)
    3672           0 :     s1 = sum(x**2)
    3673           0 :     s2 = sum(cos(c*x))
    3674           0 :     ackley = -a * exp(-b*sqrt(1.0_dp/real(n,dp)*s1)) - exp(1.0_dp/real(n,dp)*s2) + a + exp(1.0_dp)
    3675             : 
    3676           0 :   end function ackley
    3677             : 
    3678             :   ! ------------------------------------------------------------------
    3679             :   !
    3680             :   !  The Michalewicz function, N >= 2.
    3681             :   !  Search domain: x restricted to (0, Pi)
    3682             :   !  Solution:
    3683             :   !     numerical, so far best found
    3684             :   !     x(1:2)  = (/ 2.2029262967_dp, 1.5707721052_dp /)
    3685             :   !     f(x)    = -1.8013033793_dp
    3686             :   !
    3687             :   !     numerical, so far best found
    3688             :   !     x(1:5)  = (/ 2.2029054902_dp, 1.5707963436_dp, 1.2849915892_dp, 1.9230584622_dp, 1.7204697668_dp /)
    3689             :   !     f(x)    = -4.6876581791_dp
    3690             :   !     known from literature:
    3691             :   !     f(x)    = -4.687_dp
    3692             :   !
    3693             :   !     x(1:10) = (/ 2.2029055201726093_dp, 2.10505573543129_dp, &
    3694             :   !                  2.2193332517481035_dp, 1.9230584698663626_dp, &
    3695             :   !                  0.9966770382174085_dp, 2.0274797779024762_dp, &
    3696             :   !                  1.7114837946034247_dp, 1.3605717365168617_dp, &
    3697             :   !                  1.2828240065709524_dp, 1.5707963267948966_dp /)
    3698             :   !     f(x)    = -7.209069703423156_dp
    3699             :   !     known from literature:
    3700             :   !     f(x)    = -9.66_dp
    3701             :   !
    3702             :   !  Discussion:
    3703             :   !    The Michalewicz function is a multimodal test function (n! local optima).
    3704             :   !    The parameter p defines the "steepness" of the valleys or edges. Larger p leads to
    3705             :   !    more difficult search. For very large p the function behaves like a needle in the
    3706             :   !    haystack (the function values for points in the space outside the narrow peaks give
    3707             :   !    very little information on the location of the global optimum).
    3708             :   !    http://www.geatbx.com/docu/fcnindex-01.html#P150_6749
    3709             :   !
    3710             :   !    http://www.scribd.com/doc/74351406/12/Michalewicz''s-function
    3711             :   !
    3712             :   !  Author:
    3713             :   !
    3714             :   !    Matlab code by A. Hedar
    3715             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3716             :   !
    3717             :   !  Parameters:
    3718             :   !
    3719             :   !    Input, real(dp) X(N), the argument.
    3720             :   !
    3721             : 
    3722           0 :   function michalewicz(x)
    3723             : 
    3724           0 :     use mo_constants, only: pi_dp
    3725             : 
    3726             :     implicit none
    3727             : 
    3728             :     real(dp), dimension(:), intent(in) :: x
    3729             :     real(dp)                           :: michalewicz
    3730             : 
    3731             :     integer(i4) :: n
    3732             :     integer(i4) :: j
    3733           0 :     real(dp)    :: tmp
    3734             :     integer(i4), parameter :: p = 20
    3735             : 
    3736           0 :     n = size(x)
    3737           0 :     michalewicz = 0.0_dp
    3738           0 :     do j=1, n
    3739             :        ! michalewicz = michalewicz + sin(x(j)) * sin(x(j)**2 * (real(j,dp)/pi_dp))**p
    3740           0 :        tmp = x(j)*x(j) * (real(j,dp)/pi_dp)
    3741           0 :        tmp = sin(tmp)
    3742           0 :        if (abs(tmp) .lt. 1E-15) then
    3743             :           tmp = 0.0_dp
    3744             :        else
    3745           0 :           tmp = tmp**p
    3746             :        end if
    3747           0 :        tmp = sin(x(j)) * tmp
    3748           0 :        michalewicz = michalewicz + tmp
    3749             :     end do
    3750           0 :     michalewicz = -michalewicz
    3751             : 
    3752           0 :     select case(n)
    3753             :     case(2_i4)
    3754           0 :        michalewicz = michalewicz + 1.8013033793_dp
    3755             :     case(5_i4)
    3756           0 :        michalewicz = michalewicz + 4.6876581791_dp
    3757             :     end select
    3758             : 
    3759           0 :   end function michalewicz
    3760             : 
    3761             :   ! ------------------------------------------------------------------
    3762             :   !
    3763             :   !  The Booth function, N = 2.
    3764             :   !  Solution: x(1:2) = (/ 1.0_dp, 3.0_dp /)
    3765             :   !
    3766             :   !  Author:
    3767             :   !
    3768             :   !    Matlab code by A. Hedar
    3769             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3770             :   !
    3771             :   !  Parameters:
    3772             :   !
    3773             :   !    Input, real(dp) X(N), the argument.
    3774             :   !
    3775             : 
    3776           0 :   function booth(x)
    3777             : 
    3778             :     implicit none
    3779             : 
    3780             :     real(dp), dimension(:), intent(in) :: x
    3781             :     real(dp) :: booth
    3782             : 
    3783           0 :     booth = (x(1)+2.0_dp*x(2)-7.0_dp)**2 + (2.0_dp*x(1)+x(2)-5.0_dp)**2
    3784             : 
    3785           0 :   end function booth
    3786             : 
    3787             :   ! ------------------------------------------------------------------
    3788             :   !
    3789             :   !  The Hump function, N = 2.
    3790             :   !
    3791             :   !  Search Domain:
    3792             :   !     -5.0_dp <= xi <= 5.0_dp
    3793             :   !
    3794             :   !  Solution:
    3795             :   !     x(1:2) = (/ 0.08984201310031806_dp , -0.7126564030207396_dp /)     OR
    3796             :   !     x(1:2) = (/ -0.08984201310031806_dp , 0.7126564030207396_dp /)
    3797             :   !     f(x)   = 0.0_dp
    3798             :   !
    3799             :   !  Author:
    3800             :   !
    3801             :   !    Matlab code by A. Hedar
    3802             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3803             :   !
    3804             :   !  Literature:
    3805             :   !
    3806             :   !    http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO_files/Page1621.htm
    3807             :   !
    3808             :   !  Parameters:
    3809             :   !
    3810             :   !    Input, real(dp) X(N), the argument.
    3811             :   !
    3812             : 
    3813           0 :   function hump(x)
    3814             : 
    3815             :     implicit none
    3816             : 
    3817             :     real(dp), dimension(:), intent(in) :: x
    3818             :     real(dp) :: hump
    3819             : 
    3820           0 :     hump = 1.0316285_dp + 4.0_dp*x(1)**2 - 2.1_dp*x(1)**4 + x(1)**6 / 3.0_dp + x(1)*x(2) - 4.0_dp*x(2)**2 + 4.0_dp*x(2)**4
    3821             : 
    3822           0 :   end function hump
    3823             : 
    3824             :   ! ------------------------------------------------------------------
    3825             :   !
    3826             :   !  The Levy function, N >= 2.
    3827             :   !  Solution: x(1:n) = 1.0_dp
    3828             :   !
    3829             :   !  Author:
    3830             :   !
    3831             :   !    Matlab code by A. Hedar
    3832             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3833             :   !
    3834             :   !  Parameters:
    3835             :   !
    3836             :   !    Input, real(dp) X(N), the argument.
    3837             :   !
    3838             : 
    3839           0 :   function levy(x)
    3840             : 
    3841           0 :     use mo_constants, only: pi_dp
    3842             : 
    3843             :     implicit none
    3844             : 
    3845             :     real(dp), dimension(:), intent(in) :: x
    3846             :     real(dp) :: levy
    3847             : 
    3848             :     integer(i4) :: n
    3849             :     integer(i4) :: i
    3850           0 :     real(dp), dimension(size(x)) :: z
    3851             : 
    3852           0 :     n = size(x)
    3853           0 :     z = 1.0_dp+(x-1.0_dp)/4.0_dp
    3854           0 :     levy = sin(pi_dp*z(1))**2
    3855           0 :     do i=1, n-1
    3856           0 :        levy = levy + (z(i)-1.0_dp)**2 * (1.0_dp+10.0_dp*(sin(pi_dp*z(i)+1.0_dp))**2)
    3857             :     end do
    3858           0 :     levy = levy + (z(n)-1.0_dp)**2 * (1.0_dp+(sin(2.0_dp*pi_dp*z(n)))**2)
    3859             : 
    3860           0 :   end function levy
    3861             : 
    3862             :   ! ------------------------------------------------------------------
    3863             :   !
    3864             :   !  The Matyas function, N = 2.
    3865             :   !  Solution: x(1:n) = 0.0_dp
    3866             :   !
    3867             :   !  Author:
    3868             :   !
    3869             :   !    Matlab code by A. Hedar
    3870             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3871             :   !
    3872             :   !  Parameters:
    3873             :   !
    3874             :   !    Input, real(dp) X(N), the argument.
    3875             :   !
    3876             : 
    3877           0 :   function matyas(x)
    3878             : 
    3879             :     implicit none
    3880             : 
    3881             :     real(dp), dimension(:), intent(in) :: x
    3882             :     real(dp) :: matyas
    3883             : 
    3884           0 :     matyas = 0.26_dp * (x(1)**2 + x(2)**2) -0.48_dp*x(1)*x(2)
    3885             : 
    3886           0 :   end function matyas
    3887             : 
    3888             :   ! ------------------------------------------------------------------
    3889             :   !
    3890             :   !  The Perm function, N >= 4.
    3891             :   !  Solution: forall(i=1:n) x(i) = real(i,dp)
    3892             :   !
    3893             :   !  Author:
    3894             :   !
    3895             :   !    Matlab code by A. Hedar
    3896             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3897             :   !
    3898             :   !  Parameters:
    3899             :   !
    3900             :   !    Input, real(dp) X(N), the argument.
    3901             :   !
    3902             : 
    3903           0 :   function perm(x)
    3904             : 
    3905             :     implicit none
    3906             : 
    3907             :     real(dp), dimension(:), intent(in) :: x
    3908             :     real(dp) :: perm
    3909             : 
    3910             :     integer(i4) :: n
    3911             :     integer(i4) :: j, k
    3912           0 :     real(dp) :: s_in
    3913             : 
    3914           0 :     n = size(x)
    3915           0 :     perm = 0.0_dp
    3916           0 :     do k=1, n
    3917             :        s_in = 0.0_dp
    3918           0 :        do j=1, n
    3919           0 :           s_in = s_in + (real(j,dp)**k + 0.5_dp) * ((x(j)/real(j,dp))**k - 1.0_dp)
    3920             :        enddo
    3921           0 :        perm = perm + s_in**2
    3922             :     enddo
    3923             : 
    3924             : 
    3925           0 :   end function perm
    3926             : 
    3927             :   ! ------------------------------------------------------------------
    3928             :   !
    3929             :   !  The Power sum function, N = 4.
    3930             :   !  Solution:
    3931             :   !      x    = (/ 1.0000653551_dp, 2.0087089520_dp, 1.9912589253_dp, 2.9999732609_dp /)
    3932             :   !      f(x) = 0.0000000001
    3933             :   !
    3934             :   !  Author:
    3935             :   !
    3936             :   !    Matlab code by A. Hedar
    3937             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3938             :   !
    3939             :   !  Parameters:
    3940             :   !
    3941             :   !    Input, real(dp) X(N), the argument.
    3942             :   !
    3943             : 
    3944           0 :   function power_sum(x)
    3945             : 
    3946             :     implicit none
    3947             : 
    3948             :     real(dp), dimension(:), intent(in) :: x
    3949             :     real(dp) :: power_sum
    3950             : 
    3951             :     integer(i4) :: n
    3952             :     integer(i4) :: k
    3953           0 :     real(dp), dimension(4) :: b
    3954             : 
    3955           0 :     n = size(x)
    3956           0 :     b = (/ 8.0_dp, 18.0_dp, 44.0_dp, 114.0_dp /)
    3957           0 :     power_sum = 0.0_dp
    3958           0 :     do k=1, n
    3959           0 :        power_sum = power_sum + (sum(x**k) - b(k))**2
    3960             :     end do
    3961             : 
    3962           0 :   end function power_sum
    3963             : 
    3964             :   ! ------------------------------------------------------------------
    3965             :   !
    3966             :   !  Sphere model, N = 1.
    3967             :   !  Solution: x(1:n) = 0.0_dp
    3968             :   !
    3969             :   !  Author:
    3970             :   !
    3971             :   !    Matlab code by A. Hedar
    3972             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    3973             :   !
    3974             :   !  Parameters:
    3975             :   !
    3976             :   !    Input, real(dp) X(N), the argument.
    3977             :   !
    3978             : 
    3979             :   function sphere(x)
    3980             : 
    3981             :     implicit none
    3982             : 
    3983             :     real(dp), dimension(:), intent(in) :: x
    3984             :     real(dp) :: sphere
    3985             : 
    3986             :     integer(i4) :: n
    3987             :     integer(i4) :: j
    3988             : 
    3989             :     n = size(x)
    3990             :     sphere = 0.0_dp
    3991             :     do j=1, n
    3992             :        sphere = sphere + x(j)**2
    3993             :     enddo
    3994             : 
    3995           0 :   end function sphere
    3996             : 
    3997             :   ! ------------------------------------------------------------------
    3998             :   !
    3999             :   !  The sphere model
    4000             :   !  N = 2
    4001             :   !  Solution: x(1:n) = 1.0_dp
    4002             :   !
    4003             :   !  Discussion:
    4004             :   !
    4005             :   !    The function is continuous, convex, and unimodal.
    4006             :   !
    4007             :   !  Licensing:
    4008             :   !
    4009             :   !    This code is distributed under the GNU LGPL license.
    4010             :   !
    4011             :   !  Modified:
    4012             :   !
    4013             :   !    07 January 2012
    4014             :   !
    4015             :   !  Author:
    4016             :   !
    4017             :   !    John Burkardt
    4018             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4019             :   !
    4020             :   !  Reference:
    4021             :   !
    4022             :   !    Hugues Bersini, Marco Dorigo, Stefan Langerman, Gregory Seront,
    4023             :   !    Luca Gambardella,
    4024             :   !    Results of the first international contest on evolutionary optimisation,
    4025             :   !    In Proceedings of 1996 IEEE International Conference on Evolutionary
    4026             :   !    Computation,
    4027             :   !    IEEE Press, pages 611-615, 1996.
    4028             :   !
    4029             :   !  Parameters:
    4030             :   !
    4031             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4032             :   !
    4033             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4034             :   !
    4035             : 
    4036           0 :   function sphere_model_2d(x)
    4037             : 
    4038             :     implicit none
    4039             : 
    4040             :     real(dp), dimension(:,:), intent(in) :: x
    4041             :     real(dp), dimension(size(x,2)) :: sphere_model_2d
    4042             : 
    4043             :     integer(i4) :: m
    4044             :     integer(i4) :: n
    4045             :     integer(i4) :: j
    4046             : 
    4047           0 :     m = size(x,1)
    4048           0 :     n = size(x,2)
    4049           0 :     do j = 1, n
    4050           0 :        sphere_model_2d(j) = sum( ( x(1:m,j) - 1.0_dp ) ** 2 )
    4051             :     end do
    4052             : 
    4053           0 :   end function sphere_model_2d
    4054             : 
    4055             :   ! ------------------------------------------------------------------
    4056             :   !
    4057             :   !  The axis-parallel hyper-ellipsoid function
    4058             :   !  Solution: x(1:n) = 0.0_dp
    4059             :   !
    4060             :   !  Discussion:
    4061             :   !
    4062             :   !    This function is also known as the weighted sphere model.
    4063             :   !
    4064             :   !    The function is continuous, convex, and unimodal.
    4065             :   !
    4066             :   !  Licensing:
    4067             :   !
    4068             :   !    This code is distributed under the GNU LGPL license.
    4069             :   !
    4070             :   !  Modified:
    4071             :   !
    4072             :   !    18 December 2011
    4073             :   !
    4074             :   !  Author:
    4075             :   !
    4076             :   !    John Burkardt
    4077             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4078             :   !
    4079             :   !  Reference:
    4080             :   !
    4081             :   !    Marcin Molga, Czeslaw Smutnicki,
    4082             :   !    Test functions for optimization needs.
    4083             :   !
    4084             :   !  Parameters:
    4085             :   !
    4086             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4087             :   !
    4088             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4089             :   !
    4090             : 
    4091           0 :   function axis_parallel_hyper_ellips_2d(x)
    4092             : 
    4093             :     implicit none
    4094             : 
    4095             :     real(dp), dimension(:,:), intent(in) :: x
    4096             :     real(dp), dimension(size(x,2)) :: axis_parallel_hyper_ellips_2d
    4097             : 
    4098             :     integer(i4) :: m
    4099             :     integer(i4) :: n
    4100             :     integer(i4) :: j
    4101           0 :     real(dp) :: y(size(x,1))
    4102             : 
    4103             : 
    4104           0 :     m = size(x,1)
    4105           0 :     n = size(x,2)
    4106           0 :     forall(j=1:m) y(j) = real(j,dp)
    4107             : 
    4108           0 :     do j = 1, n
    4109           0 :        axis_parallel_hyper_ellips_2d(j) = sum( y(1:m) * x(1:m,j) ** 2 )
    4110             :     end do
    4111             : 
    4112           0 :   end function axis_parallel_hyper_ellips_2d
    4113             : 
    4114             :   ! ------------------------------------------------------------------
    4115             :   !
    4116             :   !  The rotated hyper-ellipsoid function
    4117             :   !  Solution: x(1:n) = 0.0_dp
    4118             :   !
    4119             :   !  Discussion:
    4120             :   !
    4121             :   !    This function is also known as the weighted sphere model.
    4122             :   !
    4123             :   !    The function is continuous, convex, and unimodal.
    4124             :   !
    4125             :   !     There is a typographical error in Molga and Smutnicki, so that the
    4126             :   !     formula for this function is given incorrectly.
    4127             :   !
    4128             :   !  Licensing:
    4129             :   !
    4130             :   !    This code is distributed under the GNU LGPL license.
    4131             :   !
    4132             :   !  Modified:
    4133             :   !
    4134             :   !    19 December 2011
    4135             :   !
    4136             :   !  Author:
    4137             :   !
    4138             :   !    John Burkardt
    4139             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4140             :   !
    4141             :   !  Reference:
    4142             :   !
    4143             :   !    Marcin Molga, Czeslaw Smutnicki,
    4144             :   !    Test functions for optimization needs.
    4145             :   !
    4146             :   !  Parameters:
    4147             :   !
    4148             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4149             :   !
    4150             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4151             :   !
    4152             : 
    4153           0 :   function rotated_hyper_ellipsoid_2d(x)
    4154             : 
    4155             :     implicit none
    4156             : 
    4157             :     real(dp), dimension(:,:), intent(in) :: x
    4158             :     real(dp), dimension(size(x,2)) :: rotated_hyper_ellipsoid_2d
    4159             : 
    4160             :     integer(i4) :: m
    4161             :     integer(i4) :: n
    4162             :     integer(i4) :: i
    4163             :     integer(i4) :: j
    4164           0 :     real(dp) :: x_sum
    4165             : 
    4166           0 :     m = size(x,1)
    4167           0 :     n = size(x,2)
    4168           0 :     do j = 1, n
    4169             : 
    4170           0 :        rotated_hyper_ellipsoid_2d(j) = 0.0_dp
    4171           0 :        x_sum = 0.0_dp
    4172             : 
    4173           0 :        do i = 1, m
    4174           0 :           x_sum = x_sum + x(i,j)
    4175           0 :           rotated_hyper_ellipsoid_2d(j) = rotated_hyper_ellipsoid_2d(j) + x_sum**2
    4176             :        end do
    4177             : 
    4178             :     end do
    4179             : 
    4180           0 :   end function rotated_hyper_ellipsoid_2d
    4181             : 
    4182             :   ! ------------------------------------------------------------------
    4183             :   !
    4184             :   !  Rosenbrock''s valley
    4185             :   !  Solution: x(1:n) = 1.0_dp
    4186             :   !
    4187             :   !  Licensing:
    4188             :   !
    4189             :   !    This code is distributed under the GNU LGPL license.
    4190             :   !
    4191             :   !  Modified:
    4192             :   !
    4193             :   !    07 January 2012
    4194             :   !
    4195             :   !  Author:
    4196             :   !
    4197             :   !    John Burkardt
    4198             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4199             :   !
    4200             :   !  Reference:
    4201             :   !
    4202             :   !    Howard Rosenbrock,
    4203             :   !    An Automatic Method for Finding the Greatest or Least Value of a Function,
    4204             :   !    Computer Journal,
    4205             :   !    Volume 3, 1960, pages 175-184.
    4206             :   !
    4207             :   !  Parameters:
    4208             :   !
    4209             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4210             :   !
    4211             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4212             :   !
    4213             : 
    4214           0 :   function rosenbrock_2d(x)
    4215             : 
    4216             :     implicit none
    4217             : 
    4218             :     real(dp), dimension(:,:), intent(in) :: x
    4219             :     real(dp), dimension(size(x,2)) :: rosenbrock_2d
    4220             : 
    4221             :     integer(i4) :: m
    4222             :     integer(i4) :: n
    4223             :     integer(i4) :: j
    4224             : 
    4225           0 :     m = size(x,1)
    4226           0 :     n = size(x,2)
    4227           0 :     do j = 1, n
    4228           0 :        rosenbrock_2d(j) = sum ( ( 1.0_dp - x(1:m,j) )**2 ) &
    4229           0 :             + sum ( ( x(2:m,j) - x(1:m-1,j) )**2 )
    4230             :     end do
    4231             : 
    4232           0 :   end function rosenbrock_2d
    4233             : 
    4234             :   ! ------------------------------------------------------------------
    4235             :   !
    4236             :   !  Rastrigin''s function
    4237             :   !  Solution: x(1:n) = 0.0_dp
    4238             :   !
    4239             :   !  Licensing:
    4240             :   !
    4241             :   !    This code is distributed under the GNU LGPL license.
    4242             :   !
    4243             :   !  Modified:
    4244             :   !
    4245             :   !    19 December 2011
    4246             :   !
    4247             :   !  Author:
    4248             :   !
    4249             :   !    John Burkardt
    4250             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4251             :   !
    4252             :   !  Reference:
    4253             :   !
    4254             :   !    Marcin Molga, Czeslaw Smutnicki,
    4255             :   !    Test functions for optimization needs.
    4256             :   !
    4257             :   !  Parameters:
    4258             :   !
    4259             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4260             :   !
    4261             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4262             :   !
    4263             : 
    4264           0 :   function rastrigin_2d(x)
    4265             : 
    4266           0 :     use mo_constants, only: pi_dp
    4267             : 
    4268             :     implicit none
    4269             : 
    4270             :     real(dp), dimension(:,:), intent(in) :: x
    4271             :     real(dp), dimension(size(x,2)) :: rastrigin_2d
    4272             : 
    4273             :     integer(i4) :: m
    4274             :     integer(i4) :: n
    4275             :     integer(i4) :: i
    4276             :     integer(i4) :: j
    4277             : 
    4278           0 :     m = size(x,1)
    4279           0 :     n = size(x,2)
    4280           0 :     do j = 1, n
    4281             : 
    4282           0 :        rastrigin_2d(j) = real( 10 * m, dp )
    4283             : 
    4284           0 :        do i = 1, m
    4285           0 :           rastrigin_2d(j) = rastrigin_2d(j) + x(i,j) ** 2 - 10.0_dp * cos( 2.0_dp * pi_dp * x(i,j) )
    4286             :        end do
    4287             : 
    4288             :     end do
    4289             : 
    4290           0 :   end function rastrigin_2d
    4291             : 
    4292             :   ! ------------------------------------------------------------------
    4293             :   !
    4294             :   !  Schwefel''s function.
    4295             :   !  Solution: x(1:n) = 420.9687_dp
    4296             :   !
    4297             :   !  Licensing:
    4298             :   !
    4299             :   !    This code is distributed under the GNU LGPL license.
    4300             :   !
    4301             :   !  Modified:
    4302             :   !
    4303             :   !    19 December 2011
    4304             :   !
    4305             :   !  Author:
    4306             :   !
    4307             :   !    John Burkardt
    4308             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4309             :   !
    4310             :   !  Reference:
    4311             :   !
    4312             :   !    Hans-Paul Schwefel,
    4313             :   !    Numerical optimization of computer models,
    4314             :   !    Wiley, 1981,
    4315             :   !    ISBN13: 978-0471099888,
    4316             :   !    LC: QA402.5.S3813.
    4317             :   !
    4318             :   !  Parameters:
    4319             :   !
    4320             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4321             :   !
    4322             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4323             :   !
    4324             : 
    4325           0 :   function schwefel_2d(x)
    4326             : 
    4327             :     implicit none
    4328             : 
    4329             :     real(dp), dimension(:,:), intent(in) :: x
    4330             :     real(dp), dimension(size(x,2)) :: schwefel_2d
    4331             : 
    4332             :     integer(i4) :: m
    4333             :     integer(i4) :: n
    4334             :     integer(i4) :: j
    4335             : 
    4336           0 :     m = size(x,1)
    4337           0 :     n = size(x,2)
    4338           0 :     do j = 1, n
    4339           0 :        schwefel_2d(j) = -sum( x(1:m,j) * sin( sqrt( abs( x(1:m,j) ) ) ) )
    4340             :     end do
    4341             : 
    4342           0 :   end function schwefel_2d
    4343             : 
    4344             :   ! ------------------------------------------------------------------
    4345             :   !
    4346             :   !  Griewank''s function
    4347             :   !  Solution: x(1:n) = 0.0_dp
    4348             :   !
    4349             :   !  Licensing:
    4350             :   !
    4351             :   !    This code is distributed under the GNU LGPL license.
    4352             :   !
    4353             :   !  Modified:
    4354             :   !
    4355             :   !    19 December 2011
    4356             :   !
    4357             :   !  Author:
    4358             :   !
    4359             :   !    John Burkardt
    4360             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4361             :   !
    4362             :   !  Reference:
    4363             :   !
    4364             :   !    Marcin Molga, Czeslaw Smutnicki,
    4365             :   !    Test functions for optimization needs.
    4366             :   !
    4367             :   !  Parameters:
    4368             :   !
    4369             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4370             :   !
    4371             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4372             :   !
    4373             : 
    4374           0 :   function griewank_2d(x)
    4375             : 
    4376             :     implicit none
    4377             : 
    4378             :     real(dp), dimension(:,:), intent(in) :: x
    4379             :     real(dp), dimension(size(x,2)) :: griewank_2d
    4380             : 
    4381             :     integer(i4) :: m
    4382             :     integer(i4) :: n
    4383             :     integer(i4) :: j
    4384           0 :     real(dp) :: y(size(x,1))
    4385             : 
    4386           0 :     m = size(x,1)
    4387           0 :     n = size(x,2)
    4388           0 :     forall(j=1:m) y(j) = real(j,dp)
    4389           0 :     y(1:m) = sqrt( y(1:m) )
    4390             : 
    4391           0 :     do j = 1, n
    4392           0 :        griewank_2d(j) = sum( x(1:m,j) ** 2 ) / 4000.0_dp &
    4393           0 :             - product( cos( x(1:m,j) / y(1:m) ) ) + 1.0_dp
    4394             :     end do
    4395             : 
    4396           0 :   end function griewank_2d
    4397             : 
    4398             :   ! ------------------------------------------------------------------
    4399             :   !
    4400             :   !  The power sum function.
    4401             :   !  Solution: x(1:n) = 0.0_dp
    4402             :   !
    4403             :   !  Licensing:
    4404             :   !
    4405             :   !    This code is distributed under the GNU LGPL license.
    4406             :   !
    4407             :   !  Modified:
    4408             :   !
    4409             :   !    19 December 2011
    4410             :   !
    4411             :   !  Author:
    4412             :   !
    4413             :   !    John Burkardt
    4414             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4415             :   !
    4416             :   !  Reference:
    4417             :   !
    4418             :   !    Marcin Molga, Czeslaw Smutnicki,
    4419             :   !    Test functions for optimization needs.
    4420             :   !
    4421             :   !  Parameters:
    4422             :   !
    4423             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4424             :   !
    4425             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4426             :   !
    4427             : 
    4428           0 :   function power_sum_2d(x)
    4429             : 
    4430             :     implicit none
    4431             : 
    4432             :     real(dp), dimension(:,:), intent(in) :: x
    4433             :     real(dp), dimension(size(x,2)) :: power_sum_2d
    4434             : 
    4435             :     integer(i4) :: m
    4436             :     integer(i4) :: n
    4437             :     integer(i4) :: j
    4438           0 :     real(dp) :: y(size(x,1))
    4439             : 
    4440           0 :     m = size(x,1)
    4441           0 :     n = size(x,2)
    4442           0 :     forall(j=1:m) y(j) = real(j,dp)
    4443           0 :     y(1:m) = y(1:m) + 1.0_dp
    4444             : 
    4445           0 :     do j = 1, n
    4446           0 :        power_sum_2d(j) = sum( abs( x(1:m,j) ) ** y(1:m) )
    4447             :     end do
    4448             : 
    4449           0 :   end function power_sum_2d
    4450             : 
    4451             :   ! ------------------------------------------------------------------
    4452             :   !
    4453             :   !  Ackley''s function
    4454             :   !  Solution: x(1:n) = 0.0_dp
    4455             :   !
    4456             :   !  Licensing:
    4457             :   !
    4458             :   !    This code is distributed under the GNU LGPL license.
    4459             :   !
    4460             :   !  Modified:
    4461             :   !
    4462             :   !    19 December 2011
    4463             :   !
    4464             :   !  Author:
    4465             :   !
    4466             :   !    John Burkardt
    4467             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4468             :   !
    4469             :   !  Reference:
    4470             :   !
    4471             :   !    Marcin Molga, Czeslaw Smutnicki,
    4472             :   !    Test functions for optimization needs.
    4473             :   !
    4474             :   !  Parameters:
    4475             :   !
    4476             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4477             :   !
    4478             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4479             :   !
    4480             : 
    4481           0 :   function ackley_2d(x)
    4482             : 
    4483           0 :     use mo_constants, only: pi_dp
    4484             : 
    4485             :     implicit none
    4486             : 
    4487             :     real(dp), dimension(:,:), intent(in) :: x
    4488             :     real(dp), dimension(size(x,2)) :: ackley_2d
    4489             : 
    4490             :     integer(i4) :: m
    4491             :     integer(i4) :: n
    4492             :     integer(i4) :: j
    4493             :     real(dp), parameter :: a = 20.0_dp
    4494             :     real(dp), parameter :: b = 0.2_dp
    4495             :     real(dp), parameter :: c = 0.2_dp
    4496             : 
    4497           0 :     m = size(x,1)
    4498           0 :     n = size(x,2)
    4499           0 :     do j = 1, n
    4500           0 :        ackley_2d(j) = -a * exp( -b * sqrt( sum( x(1:m,j)**2 ) &
    4501             :             / real( m, dp ) ) ) &
    4502             :             - exp( sum( cos( c * pi_dp * x(1:m,j) ) ) / real( m, dp ) ) &
    4503           0 :             + a + exp( 1.0_dp )
    4504             :     end do
    4505             : 
    4506           0 :   end function ackley_2d
    4507             : 
    4508             :   ! ------------------------------------------------------------------
    4509             :   !
    4510             :   !  Michalewicz''s function
    4511             :   !  Solution: x(1:n) = 0.0_dp
    4512             :   !
    4513             :   !  Licensing:
    4514             :   !
    4515             :   !    This code is distributed under the GNU LGPL license.
    4516             :   !
    4517             :   !  Modified:
    4518             :   !
    4519             :   !    19 December 2011
    4520             :   !
    4521             :   !  Author:
    4522             :   !
    4523             :   !    John Burkardt
    4524             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4525             :   !
    4526             :   !  Reference:
    4527             :   !
    4528             :   !    Marcin Molga, Czeslaw Smutnicki,
    4529             :   !    Test functions for optimization needs.
    4530             :   !
    4531             :   !  Parameters:
    4532             :   !
    4533             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4534             :   !
    4535             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4536             :   !
    4537             : 
    4538           0 :   function michalewicz_2d(x)
    4539             : 
    4540           0 :     use mo_constants, only: pi_dp
    4541             : 
    4542             :     implicit none
    4543             : 
    4544             :     real(dp), dimension(:,:), intent(in) :: x
    4545             :     real(dp), dimension(size(x,2)) :: michalewicz_2d
    4546             : 
    4547             :     integer(i4) :: m
    4548             :     integer(i4) :: n
    4549             :     integer(i4) :: j
    4550             :     integer(i4), parameter :: p = 10
    4551           0 :     real(dp) :: y(size(x,1))
    4552             : 
    4553           0 :     m = size(x,1)
    4554           0 :     n = size(x,2)
    4555           0 :     forall(j=1:m) y(j) = real(j,dp)
    4556             : 
    4557           0 :     do j = 1, n
    4558           0 :        michalewicz_2d(j) = -sum( &
    4559           0 :             sin( x(1:m,j) ) * ( sin( x(1:m,j)**2 * y(1:m) / pi_dp ) ) ** ( 2 * p ) &
    4560           0 :             )
    4561             :     end do
    4562             : 
    4563           0 :   end function michalewicz_2d
    4564             : 
    4565             :   ! ------------------------------------------------------------------
    4566             :   !
    4567             :   !  The drop wave function
    4568             :   !  Solution: x(1:n) = 0.0_dp
    4569             :   !
    4570             :   !  Licensing:
    4571             :   !
    4572             :   !    This code is distributed under the GNU LGPL license.
    4573             :   !
    4574             :   !  Modified:
    4575             :   !
    4576             :   !    07 January 2012
    4577             :   !
    4578             :   !  Author:
    4579             :   !
    4580             :   !    John Burkardt
    4581             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4582             :   !
    4583             :   !  Reference:
    4584             :   !
    4585             :   !    Marcin Molga, Czeslaw Smutnicki,
    4586             :   !    Test functions for optimization needs.
    4587             :   !
    4588             :   !  Parameters:
    4589             :   !
    4590             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4591             :   !
    4592             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4593             :   !
    4594             : 
    4595           0 :   function drop_wave_2d(x)
    4596             : 
    4597             :     implicit none
    4598             : 
    4599             :     real(dp), dimension(:,:), intent(in) :: x
    4600             :     real(dp), dimension(size(x,2)) :: drop_wave_2d
    4601             : 
    4602             :     integer(i4) :: m
    4603             :     integer(i4) :: n
    4604             :     integer(i4) :: j
    4605           0 :     real(dp) :: rsq
    4606             : 
    4607           0 :     m = size(x,1)
    4608           0 :     n = size(x,2)
    4609           0 :     do j = 1, n
    4610             : 
    4611           0 :        rsq = sum( x(1:m,j)**2 )
    4612             : 
    4613           0 :        drop_wave_2d(j) = -( 1.0_dp + cos( 12.0_dp * sqrt( rsq ) ) ) &
    4614           0 :             / ( 0.5_dp * rsq + 2.0_dp )
    4615             : 
    4616             :     end do
    4617             : 
    4618           0 :   end function drop_wave_2d
    4619             : 
    4620             :   ! ------------------------------------------------------------------
    4621             :   !
    4622             :   !  The deceptive function
    4623             :   !  Solution: forall(i=1:n) x(i) = real(i,dp)/real(n+1,dp)
    4624             :   !
    4625             :   !  Discussion:
    4626             :   !
    4627             :   !    In dimension I, the function is a piecewise linear function with
    4628             :   !    local minima at 0 and 1.0, and a global minimum at ALPHA(I) = I/(M+1).
    4629             :   !
    4630             :   !  Licensing:
    4631             :   !
    4632             :   !    This code is distributed under the GNU LGPL license.
    4633             :   !
    4634             :   !  Modified:
    4635             :   !
    4636             :   !    19 December 2011
    4637             :   !
    4638             :   !  Author:
    4639             :   !
    4640             :   !    John Burkardt
    4641             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    4642             :   !
    4643             :   !  Reference:
    4644             :   !
    4645             :   !    Marcin Molga, Czeslaw Smutnicki,
    4646             :   !    Test functions for optimization needs.
    4647             :   !
    4648             :   !  Parameters:
    4649             :   !
    4650             :   !    Input, real(dp), dimension(:,:) :: x, the arguments size (m,n), m spatial dimension, n number of arguments.
    4651             :   !
    4652             :   !    Output, real(dp), dimension(size(x,2)) :: f, the function evaluated at the arguments.
    4653             :   !
    4654             : 
    4655           0 :   function deceptive_2d(x)
    4656             : 
    4657           0 :     use mo_utils, only: le
    4658             : 
    4659             :     implicit none
    4660             : 
    4661             :     real(dp), dimension(:,:), intent(in) :: x
    4662             :     real(dp), dimension(size(x,2)) :: deceptive_2d
    4663             : 
    4664             :     integer(i4) :: m
    4665             :     integer(i4) :: n
    4666           0 :     real(dp) :: g
    4667             :     integer(i4) :: i
    4668             :     integer(i4) :: j
    4669           0 :     real(dp) :: alpha(size(x,1))
    4670             :     real(dp), parameter :: beta = 2.0_dp
    4671             :     !
    4672             :     !  I'm just choosing ALPHA in [0,1] arbitrarily.
    4673             :     !
    4674           0 :     m = size(x,1)
    4675           0 :     n = size(x,2)
    4676           0 :     do i = 1, m
    4677           0 :        alpha(i) = real( i, dp ) / real( m + 1, dp )
    4678             :     end do
    4679             : 
    4680           0 :     do j = 1, n
    4681             : 
    4682           0 :        deceptive_2d(j) = 0.0_dp
    4683             : 
    4684           0 :        do i = 1, m
    4685             : 
    4686           0 :           if ( le(x(i,j),0.0_dp) ) then
    4687           0 :              g = x(i,j)
    4688           0 :           else if ( le(x(i,j),0.8_dp * alpha(i)) ) then
    4689           0 :              g = 0.8_dp - x(i,j) / alpha(i)
    4690           0 :           else if ( le(x(i,j),alpha(i)) ) then
    4691           0 :              g = 5.0_dp * x(i,j) / alpha(i) - 4.0_dp
    4692           0 :           else if ( le(x(i,j),( 1.0_dp + 4.0_dp * alpha(i) ) / 5.0_dp) ) then
    4693           0 :              g = 1.0_dp + 5.0_dp * ( x(i,j) - alpha(i) ) / ( alpha(i) - 1.0_dp )
    4694           0 :           else if ( le(x(i,j),1.0_dp) ) then
    4695           0 :              g = 0.8_dp + ( x(i,j) - 1.0_dp ) / ( 1.0_dp - alpha(i) )
    4696             :           else
    4697           0 :              g = x(i,j) - 1.0_dp
    4698             :           end if
    4699             : 
    4700           0 :           deceptive_2d(j) = deceptive_2d(j) + g
    4701             : 
    4702             :        end do
    4703             : 
    4704           0 :        deceptive_2d(j) = deceptive_2d(j) / real( m, dp )
    4705           0 :        deceptive_2d(j) = -( deceptive_2d(j) ** beta )
    4706             : 
    4707             :     end do
    4708             : 
    4709           0 :   end function deceptive_2d
    4710             : 
    4711             :   ! ------------------------------------------------------------------
    4712             :   ! test_optimization functions of Kalyanmoy Deb
    4713             :   ! found in Deb et al. (2002), Zitzler et al. (2000) and in Matlab Central file exchange
    4714             :   ! http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/
    4715             :   ! content/TP_NSGA2/
    4716             :   ! ------------------------------------------------------------------
    4717             : 
    4718           0 :   subroutine dtlz2_3d(paraset,obj)
    4719             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    4720             : 
    4721             :     !  dimensions:
    4722             :     !      x    is usually used with 10 dimensions
    4723             :     !      f(x) is 3-dimensional
    4724             :     !  feasible domain:
    4725             :     !      x_i \in [0,1] \forall i=1,N=10
    4726             :     !  optimum:
    4727             :     !      x_i* = 0.5
    4728             :     !      x_i* \in [0,1]
    4729             :     !  comments:
    4730             :     !      problem has a spherical Pareto-optimal front
    4731             :     !      all optimal objective function values fi* must satisfy sum(fi*^2) = 1
    4732             : 
    4733             :     !  Licensing:
    4734             :     !    This original MATLAB code was covered by the following BSD License.
    4735             :     !    Copyright (c) 2011, Song Lin
    4736             :     !    All rights reserved.
    4737             :     !
    4738             :     !    This code is distributed under the GNU LGPL license.
    4739             :     !
    4740             :     !  Author:
    4741             :     !    Song Lin Jul 2011
    4742             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    4743             :     !                                  - function, dp, etc.
    4744             :     !
    4745             :     !  Reference:
    4746             :     !    Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002).
    4747             :     !        Scalable multi-objective optimization test problems (pp. 825–830).
    4748             :     !        Presented at the Congress on Evolutionary Computation (CEC 2002).
    4749             :     !        --> Eq. 9
    4750             :     !
    4751             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    4752             :     !
    4753             :     !  Parameters:
    4754             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    4755             :     !    real(dp), intent(out) :: obj       -  returns 3d objective function.
    4756             : 
    4757             :     implicit none
    4758             : 
    4759             :     real(dp), dimension(:),              intent(in)  :: paraset
    4760             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    4761             : 
    4762             :     ! local variables
    4763             :     integer(i4)                          :: ii, npara, nobj
    4764           0 :     real(dp)                             :: gm
    4765           0 :     real(dp), dimension(size(paraset,1)) :: xx
    4766           0 :     real(dp)                             :: tt
    4767             :     real(dp), parameter                  :: pi_dp = 3.141592653589793238462643383279502884197_dp
    4768             : 
    4769           0 :     npara = size(paraset)
    4770           0 :     nobj  = 3
    4771             : 
    4772           0 :     allocate(obj(nobj))
    4773             : 
    4774           0 :     obj = 0.0_dp
    4775             : 
    4776             :     gm = 0.0_dp
    4777           0 :     do ii = nobj+1, npara
    4778           0 :        gm = gm + (paraset(ii) - 0.5_dp)**2
    4779             :     end do
    4780             : 
    4781           0 :     xx = paraset * pi_dp / 2.0_dp
    4782             : 
    4783             :     ! objective 3
    4784           0 :     tt = 1.0_dp + gm
    4785           0 :     obj(nobj) = tt * sin(xx(1))
    4786             : 
    4787             :     ! objective 2
    4788           0 :     do ii = nobj-1,2,-1
    4789             :        ! print*, 'obj2: nobj-ii=',nobj-ii,'   nobj-ii+1=',nobj-ii+1
    4790           0 :        tt = tt * cos( xx(nobj-ii) )
    4791           0 :        obj(ii) = tt * sin( xx(nobj-ii+1) )
    4792             :     end do
    4793             : 
    4794             :     ! objective 1
    4795           0 :     obj(1) = tt * cos( xx(nobj-1) )
    4796             : 
    4797           0 :     where (abs(obj) < epsilon(1.0_dp))
    4798             :        obj = 0.0_dp
    4799             :     end where
    4800             : 
    4801           0 :   end subroutine dtlz2_3d
    4802             : 
    4803           0 :   subroutine dtlz2_5d(paraset,obj)
    4804             : 
    4805             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    4806             : 
    4807             :     !  dimensions:
    4808             :     !      x    is usually used with 10 dimensions
    4809             :     !      f(x) is 5-dimensional
    4810             :     !  feasible domain:
    4811             :     !      x_i \in [0,1] \forall i=1,N=10
    4812             :     !  optimum:
    4813             :     !      x_i* = 0.5
    4814             :     !      x_i* \in [0,1]
    4815             :     !  comments:
    4816             :     !      problem has a spherical Pareto-optimal front
    4817             :     !      all optimal objective function values fi* must satisfy sum(fi*^2) = 1
    4818             : 
    4819             :     !  Licensing:
    4820             :     !    This original MATLAB code was covered by the following BSD License.
    4821             :     !    Copyright (c) 2011, Song Lin
    4822             :     !    All rights reserved.
    4823             :     !
    4824             :     !    This code is distributed under the GNU LGPL license.
    4825             :     !
    4826             :     !  Author:
    4827             :     !    Song Lin Jul 2011
    4828             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    4829             :     !                                  - function, dp, etc.
    4830             :     !
    4831             :     !  Reference:
    4832             :     !    Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002).
    4833             :     !        Scalable multi-objective optimization test problems (pp. 825–830).
    4834             :     !        Presented at the Congress on Evolutionary Computation (CEC 2002).
    4835             :     !        --> Eq. 9
    4836             :     !
    4837             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    4838             :     !
    4839             :     !  Parameters:
    4840             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    4841             :     !    real(dp), intent(out) :: obj       -  returns 5d objective function.
    4842             : 
    4843             :     implicit none
    4844             : 
    4845             :     real(dp), dimension(:),              intent(in)  :: paraset
    4846             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    4847             : 
    4848             :     ! local variables
    4849             :     integer(i4)                          :: ii, npara, nobj
    4850           0 :     real(dp)                             :: gm
    4851           0 :     real(dp), dimension(size(paraset,1)) :: xx
    4852           0 :     real(dp)                             :: tt
    4853             :     real(dp), parameter                  :: pi_dp = 3.141592653589793238462643383279502884197_dp
    4854             : 
    4855           0 :     npara = size(paraset)
    4856           0 :     nobj  = 5
    4857             : 
    4858           0 :     allocate(obj(nobj))
    4859             : 
    4860           0 :     obj = 0.0_dp
    4861             : 
    4862             :     gm = 0.0_dp
    4863           0 :     do ii = nobj+1, npara
    4864           0 :        gm = gm + (paraset(ii) - 0.5)**2
    4865             :     end do
    4866             : 
    4867           0 :     xx = paraset * pi_dp / 2.0_dp
    4868             : 
    4869             :     ! objective 5
    4870           0 :     tt = 1.0_dp + gm
    4871           0 :     obj(nobj) = tt * sin(xx(1))
    4872             : 
    4873             :     ! objective 4 ... 2
    4874           0 :     do ii = nobj-1,2,-1
    4875           0 :        tt = tt * cos( xx(nobj-ii) )
    4876           0 :        obj(ii) = tt * sin( xx(nobj-ii+1) )
    4877             :     end do
    4878             : 
    4879             :     ! objective 1
    4880           0 :     obj(1) = tt * cos( xx(nobj-1) )
    4881             : 
    4882           0 :     where (abs(obj) < epsilon(1.0_dp))
    4883             :        obj = 0.0_dp
    4884             :     end where
    4885             : 
    4886           0 :   end subroutine dtlz2_5d
    4887             : 
    4888           0 :   subroutine dtlz2_10d(paraset,obj)
    4889             : 
    4890             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    4891             : 
    4892             :     !  dimensions:
    4893             :     !      x    is usually used with 10 dimensions
    4894             :     !      f(x) is 10-dimensional
    4895             :     !  feasible domain:
    4896             :     !      x_i \in [0,1] \forall i=1,N=10
    4897             :     !  optimum:
    4898             :     !      x_i* = 0.5
    4899             :     !      x_i* \in [0,1]
    4900             :     !  comments:
    4901             :     !      problem has a spherical Pareto-optimal front
    4902             :     !      all optimal objective function values fi* must satisfy sum(fi*^2) = 1
    4903             : 
    4904             :     !  Licensing:
    4905             :     !    This original MATLAB code was covered by the following BSD License.
    4906             :     !    Copyright (c) 2011, Song Lin
    4907             :     !    All rights reserved.
    4908             :     !
    4909             :     !    This code is distributed under the GNU LGPL license.
    4910             :     !
    4911             :     !  Author:
    4912             :     !    Song Lin Jul 2011
    4913             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    4914             :     !                                  - function, dp, etc.
    4915             :     !
    4916             :     !  Reference:
    4917             :     !    Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002).
    4918             :     !        Scalable multi-objective optimization test problems (pp. 825–830).
    4919             :     !        Presented at the Congress on Evolutionary Computation (CEC 2002).
    4920             :     !        --> Eq. 9
    4921             :     !
    4922             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    4923             :     !
    4924             :     !  Parameters:
    4925             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    4926             :     !    real(dp), intent(out) :: obj       -  returns 10d objective function.
    4927             : 
    4928             :     implicit none
    4929             : 
    4930             :     real(dp), dimension(:),              intent(in)  :: paraset
    4931             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    4932             : 
    4933             :     ! local variables
    4934             :     integer(i4)                          :: ii, npara, nobj
    4935           0 :     real(dp)                             :: gm
    4936           0 :     real(dp), dimension(size(paraset,1)) :: xx
    4937           0 :     real(dp)                             :: tt
    4938             :     real(dp), parameter                  :: pi_dp = 3.141592653589793238462643383279502884197_dp
    4939             : 
    4940           0 :     npara = size(paraset)
    4941           0 :     nobj  = 10
    4942             : 
    4943           0 :     allocate(obj(nobj))
    4944             : 
    4945           0 :     obj = 0.0_dp
    4946             : 
    4947             :     gm = 0.0_dp
    4948           0 :     do ii = nobj+1, npara
    4949           0 :        gm = gm + (paraset(ii) - 0.5)**2
    4950             :     end do
    4951             : 
    4952           0 :     xx = paraset * pi_dp / 2.0_dp
    4953             : 
    4954             :     ! objective 10
    4955           0 :     tt = 1.0_dp + gm
    4956           0 :     obj(nobj) = tt * sin(xx(1))
    4957             : 
    4958             :     ! objective 9 ... 2
    4959           0 :     do ii = nobj-1,2,-1
    4960           0 :        tt = tt * cos( xx(nobj-ii) )
    4961           0 :        obj(ii) = tt * sin( xx(nobj-ii+1) )
    4962             :     end do
    4963             : 
    4964             :     ! objective 1
    4965           0 :     obj(1) = tt * cos( xx(nobj-1) )
    4966             : 
    4967           0 :     where (abs(obj) < epsilon(1.0_dp))
    4968             :        obj = 0.0_dp
    4969             :     end where
    4970             : 
    4971           0 :   end subroutine dtlz2_10d
    4972             : 
    4973           0 :   subroutine fon_2d(paraset,obj)
    4974             : 
    4975             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    4976             : 
    4977             :     !  dimensions:
    4978             :     !      x    is 3 dimensional
    4979             :     !      f(x) is 2-dimensional
    4980             :     !  feasible domain:
    4981             :     !      x_i \in [-4,4] \forall i=1,N=3
    4982             :     !  optimum:
    4983             :     !      x1 = x2 = x3 \in [-1/sqrt(3), 1/sqrt(3)]
    4984             :     !  comments:
    4985             :     !      nonconvex pareto front
    4986             : 
    4987             :     !  Licensing:
    4988             :     !
    4989             :     !    This code is distributed under the GNU LGPL license.
    4990             :     !
    4991             :     !  Author:
    4992             :     !    Juliane Mai Sep 2015   - function, dp, etc.
    4993             :     !    Modified
    4994             :     !
    4995             :     !  Reference:
    4996             :     !    C. M. Fonseca and P. J. Fleming (1993)
    4997             :     !        Genetic algorithms for multiobjective optimization: Formulation, discussion and generalization
    4998             :     !        In Proceedings of the Fifth International Conference on Genetic Algorithms
    4999             :     !        S. Forrest, Ed. San Mateo, CA: Morgan Kauffman, pp. 416–423.
    5000             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5001             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5002             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5003             :     !
    5004             :     !  Parameters:
    5005             :     !    real(dp), intent(in)  :: X(:)     -  the argument of the objective function.
    5006             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5007             : 
    5008             :     implicit none
    5009             : 
    5010             :     real(dp), dimension(:),              intent(in)  :: paraset
    5011             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5012             : 
    5013             :     ! local variables
    5014             :     integer(i4) :: ii, npara, nobj
    5015           0 :     real(dp)    :: gg
    5016             : 
    5017           0 :     npara = size(paraset)
    5018           0 :     if (npara .ne. 3) stop ('mo_objective: fon_2d: This function requires 3-dimensional parameter sets')
    5019             : 
    5020           0 :     nobj  = 2
    5021           0 :     allocate(obj(nobj))
    5022             : 
    5023             :     ! objective 1
    5024           0 :     gg = 0.0_dp
    5025           0 :     do ii=1, npara
    5026           0 :        gg = gg + (paraset(ii) - 1.0/sqrt(3.0_dp))**2
    5027             :     end do
    5028           0 :     obj(1) = 1.0_dp - exp(-gg)
    5029             : 
    5030             :     ! objective 2
    5031           0 :     gg = 0.0_dp
    5032           0 :     do ii=1, npara
    5033           0 :        gg = gg + (paraset(ii) + 1.0/sqrt(3.0_dp))**2
    5034             :     end do
    5035           0 :     obj(2) = 1.0_dp - exp(-gg)
    5036             : 
    5037           0 :   end subroutine fon_2d
    5038             : 
    5039           0 :   subroutine kur_2d(paraset,obj)
    5040             : 
    5041             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5042             : 
    5043             :     !  dimensions:
    5044             :     !      x    is usually used with 3 dimensions
    5045             :     !      f(x) is 2-dimensional
    5046             :     !  feasible domain:
    5047             :     !      x_i \in [-5,5] \forall i=1,N=3
    5048             :     !  optimum:
    5049             :     !      refer to
    5050             :     !          Deb, K. (2001)
    5051             :     !              Multiobjective Optimization Using Evolutionary Algorithms.
    5052             :     !              Chichester, U.K.: Wiley.
    5053             :     !  comments:
    5054             :     !      nonconvex pareto front
    5055             : 
    5056             :     !  Licensing:
    5057             :     !    This original MATLAB code was covered by the following BSD License.
    5058             :     !    Copyright (c) 2011, Song Lin
    5059             :     !    All rights reserved.
    5060             :     !
    5061             :     !    This code is distributed under the GNU LGPL license.
    5062             :     !
    5063             :     !  Author:
    5064             :     !    Song Lin Jul 2011
    5065             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5066             :     !                                  - function, dp, etc.
    5067             :     !
    5068             :     !  Reference:
    5069             :     !    Kursawe, F. (1990)
    5070             :     !        A variant of evolution strategies for vector optimization
    5071             :     !        In Parallel Problem Solving from Nature
    5072             :     !        H.-P. Schwefel and R. Maenner, Eds.
    5073             :     !        Berlin, Germany: Springer-Verlag, pp. 193–197.
    5074             :     !    Deb, K. (2001)
    5075             :     !        Multiobjective Optimization Using Evolutionary Algorithms.
    5076             :     !        Chichester, U.K.: Wiley.
    5077             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5078             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5079             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5080             :     !
    5081             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5082             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/kur/
    5083             :     !
    5084             :     !  Parameters:
    5085             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5086             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5087             : 
    5088             :     implicit none
    5089             : 
    5090             :     real(dp), dimension(:),              intent(in)  :: paraset
    5091             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5092             : 
    5093             :     ! local variables
    5094             :     integer(i4) :: ii, npara, nobj
    5095             : 
    5096           0 :     npara = size(paraset)
    5097           0 :     nobj  = 2
    5098             : 
    5099           0 :     allocate(obj(nobj))
    5100             : 
    5101             :     ! objective 1
    5102           0 :     obj(1) = 0.0_dp
    5103           0 :     do ii=1, npara-1
    5104           0 :        obj(1) = obj(1) - 10.0_dp * exp(-0.2_dp*sqrt(paraset(ii)**2 + paraset(ii+1)**2) );
    5105             :     end do
    5106             : 
    5107             :     ! objective 2
    5108           0 :     obj(2) = 0.0_dp
    5109           0 :     do ii=1,npara
    5110           0 :        obj(2) = obj(2) + abs(paraset(ii))**0.8_dp + 5.0_dp * sin(paraset(ii)**3);
    5111             :     end do
    5112             : 
    5113           0 :   end subroutine kur_2d
    5114             : 
    5115           0 :   subroutine pol_2d(paraset,obj)
    5116             : 
    5117             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5118             : 
    5119             :     !  dimensions:
    5120             :     !      x    is 2 dimensional
    5121             :     !      f(x) is 2-dimensional
    5122             :     !  feasible domain:
    5123             :     !      x_i \in [-Pi,Pi] \forall i=1,N=2
    5124             :     !  optimum:
    5125             :     !      refer to
    5126             :     !          Deb, K. (2001)
    5127             :     !              Multiobjective Optimization Using Evolutionary Algorithms.
    5128             :     !              Chichester, U.K.: Wiley
    5129             :     !  comments:
    5130             :     !      nonconvex, disconnected pareto front
    5131             : 
    5132             :     !  Licensing:
    5133             :     !
    5134             :     !    This code is distributed under the GNU LGPL license.
    5135             :     !
    5136             :     !  Author:
    5137             :     !    Juliane Mai Sep 2015   - function, dp, etc.
    5138             :     !    Modified
    5139             :     !
    5140             :     !  Reference:
    5141             :     !    Poloni, C. (1997)
    5142             :     !        Hybrid GA for multiobjective aerodynamic shape optimization
    5143             :     !        In Genetic Algorithms in Engineering and Computer Science
    5144             :     !        G. Winter, J. Periaux, M. Galan, and P. Cuesta, Eds.
    5145             :     !        New York: Wiley, pp. 397–414.
    5146             :     !    Deb, K. (2001)
    5147             :     !        Multiobjective Optimization Using Evolutionary Algorithms.
    5148             :     !        Chichester, U.K.: Wiley
    5149             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5150             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5151             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5152             :     !
    5153             :     !  Parameters:
    5154             :     !    real(dp), intent(in)  :: X(:)     -  the argument of the objective function.
    5155             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5156             : 
    5157             :     implicit none
    5158             : 
    5159             :     real(dp), dimension(:),              intent(in)  :: paraset
    5160             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5161             : 
    5162             :     ! local variables
    5163             :     integer(i4) :: npara, nobj
    5164           0 :     real(dp)    :: a1, a2, b1, b2
    5165             : 
    5166           0 :     npara = size(paraset)
    5167           0 :     if (npara .ne. 2) stop ('mo_objective: pol_2d: This function requires 2-dimensional parameter sets')
    5168             : 
    5169           0 :     nobj  = 2
    5170           0 :     allocate(obj(nobj))
    5171             : 
    5172           0 :     a1 = 0.5_dp * sin(1.0_dp)     - 2.0_dp * cos(1.0_dp)     + 1.0_dp * sin(2.0_dp)     - 1.5_dp * cos(2.0_dp)
    5173           0 :     a2 = 1.5_dp * sin(1.0_dp)     - 1.0_dp * cos(1.0_dp)     + 2.0_dp * sin(2.0_dp)     - 0.5_dp * cos(2.0_dp)
    5174           0 :     b1 = 0.5_dp * sin(paraset(1)) - 2.0_dp * cos(paraset(1)) + 1.0_dp * sin(paraset(2)) - 1.5_dp * cos(paraset(2))
    5175           0 :     b2 = 1.5_dp * sin(paraset(1)) - 1.0_dp * cos(paraset(1)) + 2.0_dp * sin(paraset(2)) - 0.5_dp * cos(paraset(2))
    5176             : 
    5177             :     ! objective 1
    5178           0 :     obj(1) = 1.0_dp + (a1-b1)**2 + (a2-b2)**2
    5179             : 
    5180             :     ! objective 2
    5181           0 :     obj(2) = (paraset(1) + 3.0_dp)**2 + (paraset(2) + 1.0_dp)**2
    5182             : 
    5183           0 :   end subroutine pol_2d
    5184             : 
    5185           0 :   subroutine sch_2d(paraset,obj)
    5186             : 
    5187             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5188             : 
    5189             :     !  dimensions:
    5190             :     !      x    is 1-dimensional
    5191             :     !      f(x) is 2-dimensional
    5192             :     !  feasible domain:
    5193             :     !      x_i \in [-1000,1000] \forall i=1,N=1
    5194             :     !  optimum:
    5195             :     !      x \in [0,2]
    5196             :     !  comments:
    5197             :     !      convex pareto front
    5198             : 
    5199             :     !  Licensing:
    5200             :     !
    5201             :     !    This code is distributed under the GNU LGPL license.
    5202             :     !
    5203             :     !  Author:
    5204             :     !    Juliane Mai Sep 2015   - function, dp, etc.
    5205             :     !    Modified
    5206             :     !
    5207             :     !  Reference:
    5208             :     !    Schaffer J.D. (1987)
    5209             :     !        Multiple objective optimization with vector evaluated genetic algorithms
    5210             :     !        In Proceedings of the First International Conference on Genetic Algorithms,
    5211             :     !        J. J. Grefensttete, Ed. Hillsdale, NJ: Lawrence Erlbaum, pp. 93–100.
    5212             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5213             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5214             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5215             :     !
    5216             :     !  Parameters:
    5217             :     !    real(dp), intent(in)  :: X(:)     -  the argument of the objective function.
    5218             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5219             : 
    5220             :     implicit none
    5221             : 
    5222             :     real(dp), dimension(:),              intent(in)   :: paraset
    5223             :     real(dp), dimension(:), allocatable, intent(out)  :: obj
    5224             : 
    5225             :     ! local variables
    5226             :     integer(i4) :: npara, nobj
    5227             : 
    5228           0 :     npara = size(paraset)
    5229           0 :     if (npara .ne. 1) stop ('mo_objective: sch_2d: This function requires 1-dimensional parameter sets')
    5230             : 
    5231           0 :     nobj  = 2
    5232           0 :     allocate(obj(nobj))
    5233             : 
    5234             :     ! objective 1
    5235           0 :     obj(1) = paraset(1)**2
    5236             : 
    5237             :     ! objective 2
    5238           0 :     obj(2) = (paraset(1) - 2.0_dp)**2
    5239             : 
    5240           0 :   end subroutine sch_2d
    5241             : 
    5242           0 :   subroutine zdt1_2d(paraset, obj)
    5243             : 
    5244             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5245             : 
    5246             :     !  dimensions:
    5247             :     !      x    is usually used with 30 dimensions
    5248             :     !      f(x) is 2-dimensional
    5249             :     !  feasible domain:
    5250             :     !      x_i \in [0,1] \forall i=1,N=30
    5251             :     !  optimum:
    5252             :     !      x_1 \in [0,1]
    5253             :     !      x_i = 0, i=2,N=30
    5254             :     !  comments:
    5255             :     !      convex pareto front
    5256             :     !      front: f2 = 1.0 - Sqrt(f1)
    5257             : 
    5258             :     ! The schematic tradoff looks like this
    5259             :     !   /\
    5260             :     !    |
    5261             :     !  1 .
    5262             :     !    |
    5263             :     !    |
    5264             :     !    | .
    5265             :     !    |
    5266             :     !    |   .
    5267             :     !    |
    5268             :     !    |      .
    5269             :     !    |
    5270             :     !    |           .
    5271             :     !    |
    5272             :     !    |                 .
    5273             :     !    |                        .
    5274             :     !    |                                 .
    5275             :     !    |------------------------------------------.------>
    5276             :     !                                               1
    5277             : 
    5278             :     !  Licensing:
    5279             :     !    This original MATLAB code was covered by the following BSD License.
    5280             :     !    Copyright (c) 2011, Song Lin
    5281             :     !    All rights reserved.
    5282             :     !
    5283             :     !    This code is distributed under the GNU LGPL license.
    5284             :     !
    5285             :     !  Author:
    5286             :     !    Song Lin Jul 2011
    5287             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5288             :     !                                  - function, dp, etc.
    5289             :     !
    5290             :     !  Reference:
    5291             :     !    Zitzler, E., Thiele, L., & Deb, K. (2000).
    5292             :     !        Comparison of Multiobjective Evolutionary Algorithms: Empirical Results.
    5293             :     !        Evolutionary Computation, 8(2), 173–195.
    5294             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5295             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5296             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5297             :     !
    5298             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5299             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/zdt1/
    5300             :     !
    5301             :     !  Parameters:
    5302             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5303             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5304             : 
    5305             :     implicit none
    5306             : 
    5307             :     real(dp), dimension(:),              intent(in)  :: paraset
    5308             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5309             : 
    5310             :     ! local variables
    5311             :     integer(i4) :: ii, npara, nobj
    5312           0 :     real(dp)    :: gg
    5313             : 
    5314           0 :     npara = size(paraset,1)
    5315           0 :     nobj  = 2
    5316             : 
    5317           0 :     allocate(obj(nobj))
    5318             : 
    5319             :     ! objective 1
    5320           0 :     obj(1) = paraset(1)
    5321             : 
    5322             :     ! objective 2
    5323           0 :     gg = 0.0_dp
    5324           0 :     do ii = 2, npara
    5325           0 :        gg = gg + paraset(ii)
    5326             :     end do
    5327           0 :     gg = 1.0_dp + 9.0_dp * gg / real(npara-1,dp)
    5328           0 :     obj(2) = gg * (1.0_dp - sqrt(paraset(1) / gg))
    5329             : 
    5330           0 :   end subroutine zdt1_2d
    5331             : 
    5332           0 :   subroutine zdt2_2d(paraset,obj)
    5333             : 
    5334             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5335             : 
    5336             :     !  dimensions:
    5337             :     !      x    is usually used with 30 dimensions
    5338             :     !      f(x) is 2-dimensional
    5339             :     !  feasible domain:
    5340             :     !      x_i \in [0,1] \forall i=1,N=30
    5341             :     !  optimum:
    5342             :     !      x_1 \in [0,1]
    5343             :     !      x_i = 0, i=2,N=30
    5344             :     !  comments:
    5345             :     !      nonconvex pareto front
    5346             :     !      front: f2 = 1.0 - f1**2
    5347             : 
    5348             :     !  Licensing:
    5349             :     !    This original MATLAB code was covered by the following BSD License.
    5350             :     !    Copyright (c) 2011, Song Lin
    5351             :     !    All rights reserved.
    5352             :     !
    5353             :     !    This code is distributed under the GNU LGPL license.
    5354             :     !
    5355             :     !  Author:
    5356             :     !    Song Lin Jul 2011
    5357             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5358             :     !                                  - function, dp, etc.
    5359             :     !
    5360             :     !  Reference:
    5361             :     !    Zitzler, E., Thiele, L., & Deb, K. (2000).
    5362             :     !        Comparison of Multiobjective Evolutionary Algorithms: Empirical Results.
    5363             :     !        Evolutionary Computation, 8(2), 173–195.
    5364             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5365             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5366             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5367             :     !
    5368             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5369             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/zdt2/
    5370             :     !
    5371             :     !  Parameters:
    5372             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5373             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5374             : 
    5375             :     implicit none
    5376             : 
    5377             :     real(dp), dimension(:),              intent(in)  :: paraset
    5378             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5379             : 
    5380             :     ! local variables
    5381             :     integer(i4) :: npara, nobj
    5382           0 :     real(dp)    :: gg
    5383             : 
    5384           0 :     npara = size(paraset)
    5385           0 :     nobj  = 2
    5386             : 
    5387           0 :     allocate(obj(nobj))
    5388             : 
    5389             :     ! objective 1
    5390           0 :     obj(1) = paraset(1)
    5391             : 
    5392             :     ! objective 2
    5393           0 :     gg = 1.0_dp + 9.0_dp * sum(paraset(2:npara))/real(npara-1,dp)
    5394           0 :     obj(2) = gg * ( 1.0_dp - (paraset(1)/gg)**2 )
    5395             : 
    5396           0 :   end subroutine zdt2_2d
    5397             : 
    5398           0 :   subroutine zdt3_2d(paraset,obj)
    5399             : 
    5400             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5401             : 
    5402             :     !  dimensions:
    5403             :     !      x    is usually used with 30 dimensions
    5404             :     !      f(x) is 2-dimensional
    5405             :     !  feasible domain:
    5406             :     !      x_i \in [0,1] \forall i=1,N=30
    5407             :     !  optimum:
    5408             :     !      x_1 \in [0,1]
    5409             :     !      x_i = 0, i=2,N=30
    5410             :     !  comments:
    5411             :     !      convex, diconnected pareto front
    5412             :     !      front: f1 \in F and f2 = 1.0 - sqrt(f1) - f1*sin(10*Pi*f1)
    5413             :     !             where
    5414             :     !                 F = [0.0000000000, 0.0830015349] v [0.1822287280, 0.2577623634] v
    5415             :     !                     [0.4093136748, 0.4538821041] v [0.6183967944, 0.6525117038] v
    5416             :     !                     [0.8233317983, 0.8518328654]
    5417             : 
    5418             :     !  Licensing:
    5419             :     !    This original MATLAB code was covered by the following BSD License.
    5420             :     !    Copyright (c) 2011, Song Lin
    5421             :     !    All rights reserved.
    5422             :     !
    5423             :     !    This code is distributed under the GNU LGPL license.
    5424             :     !
    5425             :     !  Author:
    5426             :     !    Song Lin Jul 2011
    5427             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5428             :     !                                  - function, dp, etc.
    5429             :     !
    5430             :     !  Reference:
    5431             :     !    Zitzler, E., Thiele, L., & Deb, K. (2000).
    5432             :     !        Comparison of Multiobjective Evolutionary Algorithms: Empirical Results.
    5433             :     !        Evolutionary Computation, 8(2), 173–195.
    5434             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5435             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5436             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5437             :     !
    5438             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5439             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/zdt3/
    5440             :     !
    5441             :     !  Parameters:
    5442             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5443             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5444             : 
    5445             :     implicit none
    5446             : 
    5447             :     real(dp), dimension(:),              intent(in)  :: paraset
    5448             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5449             : 
    5450             :     ! local variables
    5451             :     integer(i4)         :: npara, nobj
    5452           0 :     real(dp)            :: gg
    5453             :     real(dp), parameter :: pi_dp = 3.141592653589793238462643383279502884197_dp
    5454             : 
    5455           0 :     npara = size(paraset)
    5456           0 :     nobj  = 2
    5457             : 
    5458           0 :     allocate(obj(nobj))
    5459             : 
    5460             :     ! objective 1
    5461           0 :     obj(1) = paraset(1)
    5462             : 
    5463             :     ! objective 2
    5464           0 :     gg = 1.0_dp + 9.0_dp * sum(paraset(2:npara))/real(npara-1,dp)
    5465           0 :     obj(2) = gg * ( 1.0_dp - sqrt(paraset(1)/gg) - paraset(1)/gg * sin(10.0_dp*pi_dp*paraset(1)) )
    5466             : 
    5467           0 :   end subroutine zdt3_2d
    5468             : 
    5469           0 :   subroutine zdt4_2d(paraset,obj)
    5470             : 
    5471             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5472             : 
    5473             :     !  dimensions:
    5474             :     !      x    is usually used with 10 dimensions
    5475             :     !      f(x) is 2-dimensional
    5476             :     !  feasible domain:
    5477             :     !      x_1 \in [0,1]
    5478             :     !      x_i \in [-5,5] \forall i=2,N=10
    5479             :     !  optimum:
    5480             :     !      x_1 \in [0,1]
    5481             :     !      x_i = 0, i=2,N=30
    5482             :     !  comments:
    5483             :     !      nonconvex pareto front
    5484             :     !      front: f2 = 1.0 - sqrt(f1)
    5485             : 
    5486             :     !  Licensing:
    5487             :     !    This original MATLAB code was covered by the following BSD License.
    5488             :     !    Copyright (c) 2011, Song Lin
    5489             :     !    All rights reserved.
    5490             :     !
    5491             :     !    This code is distributed under the GNU LGPL license.
    5492             :     !
    5493             :     !  Author:
    5494             :     !    Song Lin Jul 2011
    5495             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5496             :     !                                  - function, dp, etc.
    5497             :     !
    5498             :     !  Reference:
    5499             :     !    Zitzler, E., Thiele, L., & Deb, K. (2000).
    5500             :     !        Comparison of Multiobjective Evolutionary Algorithms: Empirical Results.
    5501             :     !        Evolutionary Computation, 8(2), 173–195.
    5502             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5503             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5504             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5505             :     !
    5506             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5507             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/zdt4/
    5508             :     !
    5509             :     !  Parameters:
    5510             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5511             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5512             : 
    5513             :     implicit none
    5514             : 
    5515             :     real(dp), dimension(:),              intent(in)  :: paraset
    5516             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5517             : 
    5518             :     ! local variables
    5519             :     integer(i4)         :: ii, npara, nobj
    5520           0 :     real(dp)            :: gg
    5521             :     real(dp), parameter :: pi_dp = 3.141592653589793238462643383279502884197_dp
    5522             : 
    5523           0 :     npara = size(paraset)
    5524           0 :     nobj  = 2
    5525             : 
    5526           0 :     allocate(obj(nobj))
    5527             : 
    5528             :     ! objective 1
    5529           0 :     obj(1) = paraset(1)
    5530             : 
    5531             :     ! objective 2
    5532           0 :     gg = 1.0 + 10.0_dp*real(npara-1,dp)
    5533           0 :     do ii = 2, npara
    5534           0 :        gg = gg + paraset(ii)**2 - 10.0_dp * cos(4.0_dp*pi_dp*paraset(ii))
    5535             :     end do
    5536           0 :     obj(2) = gg * ( 1.0_dp - sqrt(paraset(1)/gg ) )
    5537             : 
    5538           0 :   end subroutine zdt4_2d
    5539             : 
    5540           0 :   subroutine zdt6_2d(paraset,obj)
    5541             : 
    5542             :     ! Has to be a subroutine because not all compilers can handle allocatable returns of functions, e.g. gnu v4.6
    5543             : 
    5544             :     !  dimensions:
    5545             :     !      x    is usually used with 10 dimensions
    5546             :     !      f(x) is 2-dimensional
    5547             :     !  feasible domain:
    5548             :     !      x_i \in [0,1]  \forall i=1,N=10
    5549             :     !  optimum:
    5550             :     !      x_1 \in [0,1]
    5551             :     !      x_i = 0, i=2,N=30
    5552             :     !  comments:
    5553             :     !      nonconvex, nonuniformly spaced pareto front
    5554             :     !      front: f2 = 1.0 - f1**2
    5555             :     !             f1 \in [0.2807753191, 1.0]
    5556             : 
    5557             :     !  Licensing:
    5558             :     !    This original MATLAB code was covered by the following BSD License.
    5559             :     !    Copyright (c) 2011, Song Lin
    5560             :     !    All rights reserved.
    5561             :     !
    5562             :     !    This code is distributed under the GNU LGPL license.
    5563             :     !
    5564             :     !  Author:
    5565             :     !    Song Lin Jul 2011
    5566             :     !    Modified Sep 2015 Juliane Mai - translated from Matlab
    5567             :     !                                  - function, dp, etc.
    5568             :     !
    5569             :     !  Reference:
    5570             :     !    Zitzler, E., Thiele, L., & Deb, K. (2000).
    5571             :     !        Comparison of Multiobjective Evolutionary Algorithms: Empirical Results.
    5572             :     !        Evolutionary Computation, 8(2), 173–195.
    5573             :     !    Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. (2002).
    5574             :     !        A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II.
    5575             :     !        IEEE Transactions on Evolutionary Computation, 6(2), 182–197.
    5576             :     !
    5577             :     !    http://www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4/content/TP_NSGA2/
    5578             :     !    http://www.tik.ee.ethz.ch/sop/download/supplementary/testproblems/zdt6/
    5579             :     !
    5580             :     !  Parameters:
    5581             :     !    real(dp), intent(in)  :: X(:)      -  the argument of the objective function.
    5582             :     !    real(dp), intent(out) :: obj       -  returns 2d objective function.
    5583             : 
    5584             :     implicit none
    5585             : 
    5586             :     real(dp), dimension(:),              intent(in)  :: paraset
    5587             :     real(dp), dimension(:), allocatable, intent(out) :: obj
    5588             : 
    5589             :     ! local variables
    5590             :     integer(i4)         :: npara, nobj
    5591           0 :     real(dp)            :: gg
    5592             :     real(dp), parameter :: pi_dp = 3.141592653589793238462643383279502884197_dp
    5593             : 
    5594           0 :     npara = size(paraset)
    5595           0 :     nobj  = 2
    5596             : 
    5597           0 :     allocate(obj(nobj))
    5598             : 
    5599             :     ! objective 1
    5600           0 :     obj(1) = 1.0_dp - exp(-4.0_dp*paraset(1)) * sin(6.0_dp*pi_dp*paraset(1))**6
    5601             : 
    5602             :     ! objective 2
    5603           0 :     gg = 1.0_dp + 9.0_dp * (sum(paraset(2:npara))/real(npara-1,dp))**0.25_dp
    5604           0 :     obj(2) = gg * ( 1.0_dp - (obj(1) /gg)**2 )
    5605             : 
    5606           0 :   end subroutine zdt6_2d
    5607             : 
    5608             :   ! ------------------------------------------------------------------
    5609             :   !
    5610             :   !  The Ackley function, N >= 2.
    5611             :   !  Solution: x(1:n) = 0.0_dp
    5612             :   !
    5613             :   !  Author:
    5614             :   !
    5615             :   !    Matlab code by A. Hedar
    5616             :   !    Modified Jul 2012 Matthias Cuntz - function, dp, etc.
    5617             :   !
    5618             :   !  Parameters:
    5619             :   !
    5620             :   !    Input, real(dp) X(N), the argument.
    5621             :   !
    5622             : 
    5623        4585 :   function ackley_objective(parameterset, eval, arg1, arg2, arg3)
    5624             : 
    5625           0 :     use mo_constants, only: pi_dp
    5626             :     use mo_optimization_types, only : optidata_sim
    5627             : 
    5628             :     implicit none
    5629             : 
    5630             :     real(dp), intent(in), dimension(:) :: parameterset
    5631             :     procedure(eval_interface), INTENT(IN), pointer :: eval
    5632             :     real(dp), optional, intent(in) :: arg1
    5633             :     real(dp), optional, intent(out) :: arg2
    5634             :     real(dp), optional, intent(out) :: arg3
    5635             :     real(dp) :: ackley_objective
    5636             : 
    5637             :     integer(i4) :: n
    5638             :     real(dp), parameter :: a = 20.0_dp
    5639             :     real(dp), parameter :: b = 0.2_dp
    5640             :     real(dp), parameter :: c = 2.0_dp*pi_dp
    5641        4585 :     real(dp) :: s1, s2
    5642        4585 :     type(optidata_sim), dimension(:), allocatable :: et_opti
    5643             : 
    5644        4585 :     call eval(parameterset, etOptiSim=et_opti)
    5645             : 
    5646        4585 :     n = size(parameterset)
    5647      142135 :     s1 = sum(parameterset**2)
    5648      142135 :     s2 = sum(cos(c*parameterset))
    5649        4585 :     ackley_objective = -a * exp(-b*sqrt(1.0_dp/real(n,dp)*s1)) - exp(1.0_dp/real(n,dp)*s2) + a + exp(1.0_dp)
    5650             : 
    5651        4585 :   end function ackley_objective
    5652             : 
    5653      200406 :   function griewank_objective(parameterset, eval, arg1, arg2, arg3)
    5654             : 
    5655        4585 :     use mo_kind, only: i4, dp
    5656             :     use mo_optimization_types, only : optidata_sim
    5657             : 
    5658             :     implicit none
    5659             : 
    5660             :     real(dp), intent(in), dimension(:) :: parameterset
    5661             :     procedure(eval_interface), INTENT(IN), pointer :: eval
    5662             :     real(dp), optional, intent(in) :: arg1
    5663             :     real(dp), optional, intent(out) :: arg2
    5664             :     real(dp), optional, intent(out) :: arg3
    5665             :     real(dp) :: griewank_objective
    5666             : 
    5667             :     integer(i4) :: nopt
    5668             :     integer(i4) :: j
    5669      200406 :     real(dp)    :: d, u1, u2
    5670      200406 :     type(optidata_sim), dimension(:), allocatable :: et_opti
    5671             : 
    5672      200406 :     call eval(parameterset, etOptiSim=et_opti)
    5673             : 
    5674      200406 :     nopt = size(parameterset)
    5675      200406 :     if (nopt .eq. 2) then
    5676             :        d = 200.0_dp
    5677             :     else
    5678      200406 :        d = 4000.0_dp
    5679             :     end if
    5680     2201624 :     u1 = sum(parameterset**2) / d
    5681      200406 :     u2 = 1.0_dp
    5682     2201624 :     do j=1, nopt
    5683     2201624 :        u2 = u2 * cos(parameterset(j)/sqrt(real(j,dp)))
    5684             :     end do
    5685      200406 :     griewank_objective = u1 - u2 + 1.0_dp
    5686             :     !
    5687      200406 :   end function griewank_objective
    5688             : 
    5689             : 
    5690     1492502 :   subroutine eval_dummy(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim, &
    5691             :     lake_level, lake_volume, lake_area, lake_spill, lake_outflow, BFI)
    5692      200406 :     use mo_kind, only : dp
    5693             :     use mo_optimization_types, only : optidata_sim, optidata
    5694             : 
    5695             :     implicit none
    5696             : 
    5697             :     real(dp),    dimension(:), intent(in) :: parameterset
    5698             :     integer(i4), dimension(:),                 optional, intent(in)  :: opti_domain_indices
    5699             :     real(dp),    dimension(:, :), allocatable, optional, intent(out) :: runoff        ! dim1=time dim2=gauge
    5700             :     type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim       ! dim1=ncells, dim2=time
    5701             :     type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim ! dim1=ncells, dim2=time
    5702             :     type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim       ! dim1=ncells, dim2=time
    5703             :     type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim      ! dim1=ncells, dim2=time
    5704             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_level    !< dim1=time dim2=lake
    5705             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_volume   !< dim1=time dim2=lake
    5706             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_area     !< dim1=time dim2=lake
    5707             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_spill    !< dim1=time dim2=lake
    5708             :     real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_outflow  !< dim1=time dim2=lake
    5709             :   real(dp),    dimension(:), allocatable, optional, intent(out) :: BFI         !< baseflow index, dim1=domainID
    5710             : 
    5711     1492502 :     type(optidata) :: dummyData
    5712             :     integer(i4) :: i
    5713             : 
    5714     1492502 :     allocate(dummyData%dataObs(1, 1))
    5715     5970008 :     dummyData%dataObs = 0.0_dp
    5716             : 
    5717     1492502 :     if (present(etOptiSim)) then
    5718           0 :       do i=1, size(etOptiSim)
    5719           0 :         call etOptiSim(i)%init(dummyData)
    5720             :       end do
    5721             :     end if
    5722             : 
    5723     1492502 :     if (present(neutronsOptiSim)) then
    5724           0 :       do i=1, size(neutronsOptiSim)
    5725           0 :         call neutronsOptiSim(1)%init(dummyData)
    5726             :       end do
    5727             :     end if
    5728             : 
    5729     1492502 :     if (present(twsOptiSim)) then
    5730           0 :       do i=1, size(twsOptiSim)
    5731           0 :         call twsOptiSim(1)%init(dummyData)
    5732             :       end do
    5733             :     end if
    5734             : 
    5735     1492502 :     if (present(smOptiSim)) then
    5736           0 :       do i=1, size(smOptiSim)
    5737           0 :         call smOptiSim(1)%init(dummyData)
    5738             :       end do
    5739             :     end if
    5740             : 
    5741     1492502 :     if (present(runoff)) then
    5742           0 :       allocate(runoff(1, 1))
    5743           0 :       runoff(:, :) = 0.0_dp
    5744             :     end if
    5745             : 
    5746     1492502 :     if (present(BFI)) then
    5747           0 :       allocate(BFI(1))
    5748           0 :       BFI(:) = 0.0_dp
    5749             :     end if
    5750             : 
    5751     1492502 :     if (present(lake_level)) then
    5752           0 :       allocate(lake_level(1, 1))
    5753           0 :       lake_level(:, :) = 0.0_dp
    5754             :     end if
    5755             : 
    5756     1492502 :     if (present(lake_volume)) then
    5757           0 :       allocate(lake_volume(1, 1))
    5758           0 :       lake_volume(:, :) = 0.0_dp
    5759             :     end if
    5760             : 
    5761     1492502 :     if (present(lake_area)) then
    5762           0 :       allocate(lake_area(1, 1))
    5763           0 :       lake_area(:, :) = 0.0_dp
    5764             :     end if
    5765             : 
    5766     1492502 :     if (present(lake_spill)) then
    5767           0 :       allocate(lake_spill(1, 1))
    5768           0 :       lake_spill(:, :) = 0.0_dp
    5769             :     end if
    5770             : 
    5771     1492502 :     if (present(lake_outflow)) then
    5772           0 :       allocate(lake_outflow(1, 1))
    5773           0 :       lake_outflow(:, :) = 0.0_dp
    5774             :     end if
    5775             : 
    5776     1492502 :   end subroutine eval_dummy
    5777             : 
    5778             : END MODULE mo_opt_functions

Generated by: LCOV version 1.16