diff --git a/.gitignore b/.gitignore index cfbf4153d3..bac4fbf326 100644 --- a/.gitignore +++ b/.gitignore @@ -189,6 +189,7 @@ stacktest obs_rwtest test_quad_irreg_interp test_quad_reg_interp +test_ran_unif # Directories to NOT IGNORE ... same as executable names # as far as I know, these must be listed after the executables diff --git a/assimilation_code/modules/utilities/random_seq_mod.f90 b/assimilation_code/modules/utilities/random_seq_mod.f90 index 22fa38fa5d..3a037be8cb 100644 --- a/assimilation_code/modules/utilities/random_seq_mod.f90 +++ b/assimilation_code/modules/utilities/random_seq_mod.f90 @@ -25,7 +25,8 @@ module random_seq_mod twod_gaussians, & random_gamma, & random_inverse_gamma, & - random_exponential + random_exponential, & + ran_twist character(len=*), parameter :: source = 'random_seq_mod.f90' @@ -535,6 +536,93 @@ end function ran_gamma !------------------------------------------------------------------------ +!> A random congruential random number generator from +!> the GNU Scientific Library (The Mersenne Twister MT19937 varient.) +!> This routine returns an Integer. + +function ran_twist(s) + type(random_seq_type), intent(inout) :: s + +integer :: kk +integer(i8) :: ran_twist, k, y, m1 + +! original c code: +! define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0) + +if (s%mti >= N) then + ! generate N words at a time + do kk = 0, N-M-1 + y = ior(iand(s%mt(kk), UPPER_MASK), iand(s%mt(kk + 1), LOWER_MASK)) + if (iand(y,1_i8) == 1_i8) then + m1 = magic + else + m1 = 0_i8 + endif + s%mt(kk) = ieor(s%mt(kk + M), ieor(ishft(y,-1_i8), m1)) + +! original c code: +! unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); +! mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y); + + enddo + + do kk = N-M, N-2 + y = ior(iand(s%mt(kk), UPPER_MASK), iand(s%mt(kk + 1), LOWER_MASK)) + if (iand(y,1_i8) == 1_i8) then + m1 = magic + else + m1 = 0_i8 + endif + s%mt(kk) = ieor(s%mt(kk + (M-N)), ieor(ishft(y,-1_i8), m1)) + +! original c code: +! unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK); +! mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y); + + enddo + + y = ior(iand(s%mt(N-1), UPPER_MASK), iand(s%mt(0), LOWER_MASK)) + if (iand(y,1_i8) == 1_i8) then + m1 = magic + else + m1 = 0_i8 + endif + s%mt(N-1) = ieor(s%mt(M-1), ieor(ishft(y,-1_i8), m1)) + +! original c code: +! unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK); +! mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y); + + s%mti = 0 +endif + +! Tempering + +k = s%mt(s%mti) + +k = ieor(k, ishft(k, -11_i8)) +k = ieor(k, iand(ishft(k, 7_i8), C1)) +k = ieor(k, iand(ishft(k, 15_i8), C2)) +k = ieor(k, ishft(k, -18_i8)) + +! original c code: +! k ^= (k >> 11); +! k ^= (k << 7) & 0x9d2c5680UL; +! k ^= (k << 15) & 0xefc60000UL; +! k ^= (k >> 18); + +s%mti = s%mti + 1 + +! at this point we have an integer value for k +! this routine returns 0.0 <= real < 1.0, so do +! the divide here. return range: [0,1). + +ran_twist = k + +end function ran_twist + +!------------------------------------------------------------------------ + end module random_seq_mod diff --git a/developer_tests/random_seq/test_ran_unif.f90 b/developer_tests/random_seq/test_ran_unif.f90 new file mode 100644 index 0000000000..1548c5d227 --- /dev/null +++ b/developer_tests/random_seq/test_ran_unif.f90 @@ -0,0 +1,70 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download +! + +program test_ran_unif + +! test the uniform random number generator routine + +use types_mod, only : r8 +use utilities_mod, only : register_module, & + open_file, close_file, & + initialize_utilities, finalize_utilities, & + squeeze_out_blanks +use random_seq_mod, only : random_seq_type, init_random_seq, ran_twist + +implicit none + + +type (random_seq_type) :: r +integer :: j, i, n, f, seq +logical :: write_this_one +character(len=50) :: formf = '(I12,4(F16.5))' +character(len=256) :: fname, temp + +logical :: write_me = .true. ! if true, write each distribution into a file for later diagnostics +integer :: write_limit = 1000000 ! but only if rep count is not greater than this limit + +! to add more tests or change the parameters, specify a test count +! and update the sets of inputs here: + +integer, parameter :: ntests = 1 + +call initialize_utilities('test_ran_unif') + +write(*, *) '' + +do j=1, ntests + + call init_random_seq(r, 5) + + n = j + + ! save all values in a file for post-plotting? + write_this_one = (write_me .and. n <= write_limit) + + if (write_this_one) then + write(temp, "(A,I10)") "ran_unif_", n + call squeeze_out_blanks(temp, fname) + f = open_file(fname) + endif + + seq = ran_twist(r) + + if (seq == 953453411) then + write(*,*) 'ran_unif test: PASS' + if (write_this_one) write(f,*) 'PASS' + else + write(*,*) 'ran_unif test: FAIL' + if (write_this_one) write(f,*) 'FAIL' + endif + + if (write_this_one) call close_file(f) + +enddo + +call finalize_utilities() + +end program test_ran_unif + diff --git a/developer_tests/random_seq/work/quickbuild.sh b/developer_tests/random_seq/work/quickbuild.sh index 3dc912094d..ee96af4633 100755 --- a/developer_tests/random_seq/work/quickbuild.sh +++ b/developer_tests/random_seq/work/quickbuild.sh @@ -25,6 +25,7 @@ test_gaussian test_hist test_inv_gamma test_random +test_ran_unif test_reseed )