diff --git a/README.md b/README.md index d341e99b0e92bbc139e5ab86d7c850af17fc8b23..52e0b10892499a6dc4c9e3517e00638acaa2c219 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,6 @@ The `libfortran-support` library includes some general Fortran supporting module - `mo_io_units`: io unit definitions - `mo_namelist`: open/close namelist files - `mo_octree`: octree data structure and operations -- `mo_random_number_generators`: generators for pseudo random numbers (should be moved to math-support once available) - `mo_simple_dump`: array value dumping - `mo_util_backtrace`: function backtrace - `mo_util_file`: file operations diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 18d805455b17878bbd4fccf5e0e5f4ef4d6061ab..fc3c7d64c4c88f0a46867e8ed20a2d9127c8f0bf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -34,7 +34,6 @@ add_library( mo_io_units.F90 mo_namelist.F90 mo_octree.F90 - mo_random_number_generators.F90 mo_simple_dump.F90 mo_util_backtrace.F90 mo_util_file.F90 diff --git a/src/mo_random_number_generators.F90 b/src/mo_random_number_generators.F90 deleted file mode 100644 index e77e04acbde284a696883202b50a4af8d9df76e0..0000000000000000000000000000000000000000 --- a/src/mo_random_number_generators.F90 +++ /dev/null @@ -1,152 +0,0 @@ -! ICON -! -! --------------------------------------------------------------- -! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss -! Contact information: icon-model.org -! -! See AUTHORS.TXT for a list of authors -! See LICENSES/ for license information -! SPDX-License-Identifier: BSD-3-Clause -! --------------------------------------------------------------- - -MODULE mo_random_number_generators - - USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64 - USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32, int64, real64 - - IMPLICIT NONE - PRIVATE - - PUBLIC :: random_normal_values, clip, initialize_random_number_generator, generate_uniform_random_number - - INTERFACE generate_uniform_random_number - MODULE PROCEDURE generate_uniform_random_number_wp - END INTERFACE generate_uniform_random_number - - INTERFACE random_normal_values - MODULE PROCEDURE random_normal_values_wp - END INTERFACE random_normal_values - - INTERFACE clip - MODULE PROCEDURE clip_wp - END INTERFACE clip - - ! integer parameters for using revised minstd parameters as in - ! https://en.wikipedia.org/wiki/Lehmer_random_number_generator#Parameters_in_common_use - INTEGER(KIND=int64), PARAMETER :: minstd_A0 = 16807 - INTEGER(KIND=int64), PARAMETER :: minstd_A = 48271 - - INTEGER(KIND=int64), PARAMETER :: minstd_M = 2147483647 - REAL(KIND=real64), PARAMETER :: minstd_Scale = 1.0_real64/REAL(minstd_M, real64) - -CONTAINS - - ! initializes Lehmer random number generator - PURE FUNCTION initialize_random_number_generator(iseed, jseed) RESULT(minstd_state) - - INTEGER(KIND=int32), INTENT(IN) :: iseed - INTEGER(KIND=int32), INTENT(IN) :: jseed - - INTEGER(KIND=int32) :: minstd_state - - ! doing kind of an xorshift32 see https://en.wikipedia.org/wiki/Xorshift - minstd_state = iseed - minstd_state = IEOR(minstd_state, jseed*2**13) - minstd_state = IEOR(minstd_state, jseed/2**17) - minstd_state = IEOR(minstd_state, jseed*2**5) - - minstd_state = MOD(minstd_A0*minstd_state, minstd_M) - - minstd_state = MOD(minstd_A*minstd_state, minstd_M) - - END FUNCTION initialize_random_number_generator - - ! generated Lehmer random number generator https://en.wikipedia.org/wiki/Lehmer_random_number_generator - FUNCTION generate_uniform_random_number_wp(minstd_state) RESULT(random_number) - - INTEGER(KIND=int32), INTENT(INOUT) :: minstd_state - REAL(KIND=real64) :: random_number - - minstd_state = MOD(minstd_A*minstd_state, minstd_M) - - random_number = ABS(minstd_Scale*minstd_state) - - END FUNCTION generate_uniform_random_number_wp - - SUBROUTINE clip_wp(values_range, values) - - ! Subroutine arguments (in/out/inout) - REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers - REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers - - ! Limit random number range - WHERE (values(:) > values_range) - values(:) = values_range - ELSEWHERE(values(:) < -values_range) - values(:) = -values_range - END WHERE - - END SUBROUTINE clip_wp - - SUBROUTINE random_normal_values_wp(seed, values_range, values, mean, stdv) - - ! Subroutine arguments (in/out/inout) - INTEGER, INTENT(IN) :: seed ! seed to be used - REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers - REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers - REAL(wp), INTENT(IN), OPTIONAL :: mean - REAL(wp), INTENT(IN), OPTIONAL :: stdv - - REAL(wp), PARAMETER :: pi = 3.141592653589793238_wp ! replace with math constants from math-support once available - - ! Local variables - INTEGER(KIND=int32) :: i, rng_state, n - REAL(KIND=wp) :: u1, u2, z0, z1, magnitude, theta, mu, sigma - - IF (PRESENT(mean)) THEN - mu = mean - ELSE - mu = 0.0_wp - END IF - - IF (PRESENT(stdv)) THEN - sigma = stdv - ELSE - sigma = 1.0_wp - END IF - - n = SIZE(values) - - ! ----------------------------------------------------------------------- - ! Initialize a random number generator, by using the MINSTD linear - ! congruential generator (LCG). "nmaxstreams" indicates that random - ! numbers will be requested in blocks of this length. The generator - ! is seeded with "seed". - !------------------------------------------------------------------------- - - DO i = 1, n/2 + 1 - rng_state = initialize_random_number_generator(seed, i) - - u1 = generate_uniform_random_number(rng_state) - DO WHILE (u1 <= EPSILON(0.0_wp)) - u1 = generate_uniform_random_number(rng_state) - END DO - - u2 = generate_uniform_random_number(rng_state) - - ! using the Box-Muller transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform - magnitude = sigma*SQRT(-2.0_wp*LOG(u1)); - theta = 2.0_wp*pi*u2 - - z0 = magnitude*COS(theta) + mu; - z1 = magnitude*SIN(theta) + mu; - values(2*i - 1) = z0 - IF (2*i <= n) values(2*i) = z1 - - END DO - - CALL clip(values_range, values) - - END SUBROUTINE random_normal_values_wp - -END MODULE mo_random_number_generators diff --git a/test/fortran/helpers.f90 b/test/fortran/helpers.f90 index 602f26231fe7a776694d1c84a6f26f88e4a58a29..6c260baa2fdab096591be214f3407fcafa26ea99 100644 --- a/test/fortran/helpers.f90 +++ b/test/fortran/helpers.f90 @@ -53,56 +53,4 @@ CONTAINS END SUBROUTINE open_logfile - FUNCTION calculate_mean_wp(array) RESULT(mean) - REAL(wp), INTENT(IN) :: array(:) - REAL(wp) :: mean - INTEGER :: i, n - - mean = 0.0_wp - variance = 0.0_wp - - n = SIZE(array) - - DO i = 1, SIZE(array) - mean = mean + array(i) - END DO - mean = mean/n - - END FUNCTION - - FUNCTION calculate_variance_wp(array) RESULT(variance) - REAL(wp), INTENT(IN) :: array(:) - REAL(wp) :: mean, variance - INTEGER :: i, n - - mean = 0.0_wp - variance = 0.0_wp - - n = SIZE(array) - mean = calculate_mean_wp(array) - - DO i = 1, SIZE(array) - variance = variance + (array(i) - mean)**2 - END DO - variance = variance/n - - END FUNCTION - - SUBROUTINE assert_statistics(test_name, array, mean, variance, tol, max, min) - USE FORTUTF - CHARACTER(LEN=*), INTENT(IN) :: test_name - REAL(wp), INTENT(IN) :: array(:) - REAL(wp), INTENT(IN) :: mean, variance, tol, max, min - - CALL TAG_TEST(test_name//"__mean") - CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), mean, tol) - CALL TAG_TEST(test_name//"__variance") - CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), variance, tol) - CALL TAG_TEST(test_name//"__max") - CALL ASSERT_LESS_THAN_EQUAL(MAXVAL(array), max) - CALL TAG_TEST(test_name//"__min") - CALL ASSERT_GREATER_THAN_EQUAL(MINVAL(array), min) - - END SUBROUTINE - END MODULE helpers diff --git a/test/fortran/test_random_number_generators.f90 b/test/fortran/test_random_number_generators.f90 deleted file mode 100644 index 78f4932561920490df80d85f4390f6ad34c4e636..0000000000000000000000000000000000000000 --- a/test/fortran/test_random_number_generators.f90 +++ /dev/null @@ -1,101 +0,0 @@ -! ICON -! -! --------------------------------------------------------------- -! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss -! Contact information: icon-model.org -! -! See AUTHORS.TXT for a list of authors -! See LICENSES/ for license information -! SPDX-License-Identifier: BSD-3-Clause -! --------------------------------------------------------------- - -MODULE TEST_mo_random_number_generators - USE FORTUTF - USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64 - USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32 - USE helpers, ONLY: assert_statistics, calculate_mean_wp, calculate_variance_wp - -CONTAINS - SUBROUTINE TEST_CALCULATE_STATISTICS_WP - REAL(wp) :: array(8) = (/2.0_wp, 4.0_wp, 4.0_wp, 4.0_wp, 5.0_wp, 5.0_wp, 7.0_wp, 9.0_wp/) - CALL TAG_TEST("TEST_calculate_mean_wp") - CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), 5.0_wp) - CALL TAG_TEST("TEST_calculate_variance_wp") - CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), 4.0_wp) - END SUBROUTINE - - SUBROUTINE TEST_CLIP_WP - USE mo_random_number_generators - REAL(wp) :: array(6) = (/4.1_wp, -4.1_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/) - CALL TAG_TEST("TEST_clip_wp") - CALL clip(2.0_wp, array) - CALL ASSERT_EQUAL(array, (/2.0_wp, -2.0_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/)) - END SUBROUTINE - - SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC - USE mo_random_number_generators - INTEGER, PARAMETER :: n = 1e5 - INTEGER :: i - INTEGER(int32) :: rng_state - REAL(wp) :: array(n) - - rng_state = initialize_random_number_generator(11, 17) - DO i = 1, n - array(i) = generate_uniform_random_number(rng_state) - END DO - - CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC", & - array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp) - - END SUBROUTINE - - SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED - USE mo_random_number_generators - INTEGER, PARAMETER :: n = 1e5 - INTEGER :: i - INTEGER(int32) :: rng_state - REAL(wp) :: array(n) - - DO i = 1, n - rng_state = initialize_random_number_generator(i, 11) - array(i) = generate_uniform_random_number(rng_state) - END DO - - CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED", & - array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp) - - END SUBROUTINE - - SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED - USE mo_random_number_generators - INTEGER, PARAMETER :: n = 1e5 - INTEGER :: i - INTEGER(int32) :: rng_state - REAL(wp) :: array(n) - - DO i = 1, n - rng_state = initialize_random_number_generator(11, i) - array(i) = generate_uniform_random_number(rng_state) - END DO - - CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED", & - array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp) - - END SUBROUTINE - - SUBROUTINE TEST_RANDOM_NORMAL_VALUES_WP - USE mo_random_number_generators - INTEGER, PARAMETER :: n = 1e5 + 1 - REAL(wp) :: array(n) - REAL(wp) :: mean = -5.0_wp - REAL(wp) :: variance = 7.1_wp - REAL(wp) :: value_range = 1.0e32_wp - - CALL random_normal_values(seed=11, values_range=value_range, values=array, mean=mean, stdv=SQRT(variance)) - - CALL assert_statistics("TEST_RANDOM_NORMAL_VALUES_WP", & - array, mean=mean, variance=variance, tol=1.0e-2_wp, max=value_range, min=-value_range) - - END SUBROUTINE - -END MODULE TEST_mo_random_number_generators