ran_unif Mersenne twister test: Resolution to Issue #499

.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

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
 

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
+

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
 )