[Dart-dev] [5609] DART/branches/development/random_seq: Revised random number generator based on the Mersenne Twister algorithm
nancy at ucar.edu
nancy at ucar.edu
Mon Mar 19 08:50:48 MDT 2012
Revision: 5609
Author: nancy
Date: 2012-03-19 08:50:48 -0600 (Mon, 19 Mar 2012)
Log Message:
-----------
Revised random number generator based on the Mersenne Twister algorithm
from the GNU Scientific Library. It seems to have good behavior if
reseeded frequently, which is a possible usage pattern if perfect_model_obs
is called each timestep instead of running continuously. the files in
this directory incorporate the old contents of the random_nr files, so
it can be removed and all references to it in the path_names_* files can
be removed.
Modified Paths:
--------------
DART/branches/development/random_seq/random_seq_mod.f90
DART/branches/development/random_seq/random_seq_mod.html
Added Paths:
-----------
DART/branches/development/random_seq/test/
DART/branches/development/random_seq/test/input.nml
DART/branches/development/random_seq/test/mkmf_test_corr
DART/branches/development/random_seq/test/mkmf_test_diff
DART/branches/development/random_seq/test/mkmf_test_random_gsl
DART/branches/development/random_seq/test/mkmf_test_reseed
DART/branches/development/random_seq/test/path_names_test_corr
DART/branches/development/random_seq/test/path_names_test_diff
DART/branches/development/random_seq/test/path_names_test_random_gsl
DART/branches/development/random_seq/test/path_names_test_reseed
DART/branches/development/random_seq/test_random_gsl.f90
DART/branches/development/random_seq/test_reseed.f90
Removed Paths:
-------------
DART/branches/development/random_seq/input.nml
DART/branches/development/random_seq/mkmf_test_diff
DART/branches/development/random_seq/path_names_test_diff
-------------- next part --------------
Deleted: DART/branches/development/random_seq/input.nml
===================================================================
--- DART/branches/development/random_seq/input.nml 2012-03-16 23:11:46 UTC (rev 5608)
+++ DART/branches/development/random_seq/input.nml 2012-03-19 14:50:48 UTC (rev 5609)
@@ -1,10 +0,0 @@
-&utilities_nml
- TERMLEVEL = 1,
- logfilename = 'dart_log.out' /
-
-&test_diff_nml
- mean = 0.0d0,
- sd1 = 1.0d0,
- sd2 = 1.0d0,
- iseed = 123 /
-
Deleted: DART/branches/development/random_seq/mkmf_test_diff
===================================================================
--- DART/branches/development/random_seq/mkmf_test_diff 2012-03-16 23:11:46 UTC (rev 5608)
+++ DART/branches/development/random_seq/mkmf_test_diff 2012-03-19 14:50:48 UTC (rev 5609)
@@ -1,18 +0,0 @@
-#!/bin/csh
-#
-# DART software - Copyright 2004 - 2011 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
-#
-# $Id$
-
-../mkmf/mkmf -p test_diff -t ../mkmf/mkmf.template -c"-Duse_netCDF" \
- -a "../" path_names_test_diff
-
-exit $status
-
-# <next few lines under version control, do not edit>
-# $URL$
-# $Revision$
-# $Date$
-
Deleted: DART/branches/development/random_seq/path_names_test_diff
===================================================================
--- DART/branches/development/random_seq/path_names_test_diff 2012-03-16 23:11:46 UTC (rev 5608)
+++ DART/branches/development/random_seq/path_names_test_diff 2012-03-19 14:50:48 UTC (rev 5609)
@@ -1,7 +0,0 @@
-random_seq/test_diff.f90
-random_seq/random_seq_mod.f90
-random_nr/random_nr_mod.f90
-common/types_mod.f90
-utilities/utilities_mod.f90
-mpi_utilities/null_mpi_utilities_mod.f90
-time_manager/time_manager_mod.f90
Modified: DART/branches/development/random_seq/random_seq_mod.f90
===================================================================
--- DART/branches/development/random_seq/random_seq_mod.f90 2012-03-16 23:11:46 UTC (rev 5608)
+++ DART/branches/development/random_seq/random_seq_mod.f90 2012-03-19 14:50:48 UTC (rev 5609)
@@ -2,6 +2,9 @@
! provided by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
+! This module now contains both the original contents and the routines
+! which used to be in a separate random_nr module.
+
module random_seq_mod
! <next few lines under version control, do not edit>
@@ -10,9 +13,8 @@
! $Revision$
! $Date$
-use types_mod, only : r8
-use utilities_mod, only : register_module
-use random_nr_mod, only : random_seq_type, init_ran1, ran1, gasdev
+use types_mod, only : r8, digits12, i8
+use utilities_mod, only : register_module, error_handler, E_ERR
implicit none
private
@@ -27,7 +29,7 @@
revdate = "$Date$"
! Gives ability to generate unique repeatable sequences of random numbers
-! using random congruential package. Needed to allow different assim algorithms
+! using random congruential package. Needed to allow different assim algorithms
! that require random numbers to see identical observational sequences.
! Used to give different sequences a different but repeatable start
@@ -36,17 +38,36 @@
integer :: seq_number = -1
+! the following routines were transcribed from C to F90, originally
+! from the GNU scientific library: init_ran, ran_unif, ran_gauss
+
+integer, parameter :: N = 624 ! period parameters
+integer, parameter :: M = 397
+
+! hexadecimal constants
+integer(i8), parameter :: UPPER_MASK = z'80000000'
+integer(i8), parameter :: LOWER_MASK = z'7FFFFFFF'
+integer(i8), parameter :: magic = z'9908b0df'
+
+type random_seq_type
+ private
+ integer :: mti
+ integer(i8) :: mt(0:N-1) ! FIXME: move this to 1:N once its working
+ real(r8) :: lastg
+ logical :: gset
+end type random_seq_type
+
logical, save :: module_initialized = .false.
+
contains
+
!========================================================================
subroutine init_random_seq(r, seed)
!----------------------------------------------------------------------
-! You cannot generate any random numbers without calling this,
-! so this is a sufficient entry point for initializing the module.
! An integer seed can be used to get a particular repeatable sequence.
!
@@ -55,16 +76,11 @@
type(random_seq_type), intent(inout) :: r
integer, optional, intent(in) :: seed
-if ( .not. module_initialized ) then
- call register_module(source, revision, revdate)
- module_initialized = .true.
-endif
-
! Initialize the generator; use given seed if present, else sequence
if(present(seed)) then
- call init_ran1(r, seed)
+ call init_ran(r, seed)
else
- call init_ran1(r, seq_number)
+ call init_ran(r, seq_number)
seq_number = seq_number - 1
endif
@@ -79,7 +95,7 @@
type(random_seq_type), intent(inout) :: r
real(r8) :: random_uniform
-random_uniform = ran1(r)
+random_uniform = ran_unif(r)
end function random_uniform
@@ -93,7 +109,7 @@
real(r8), intent(in) :: mean, standard_deviation
real(r8) :: random_gaussian
-random_gaussian = gasdev(r) * standard_deviation + mean
+random_gaussian = ran_gauss(r) * standard_deviation + mean
end function random_gaussian
@@ -111,7 +127,7 @@
integer :: i
do i = 1, n
- rnum(i) = gasdev(r) * standard_deviation + mean
+ rnum(i) = ran_gauss(r) * standard_deviation + mean
end do
end subroutine several_random_gaussians
@@ -136,8 +152,8 @@
a22 = sqrt(cov(2, 2) - a21**2)
! Two base independent gaussian deviates
-x1 = gasdev(r)
-x2 = gasdev(r)
+x1 = ran_gauss(r)
+x2 = ran_gauss(r)
! Use these to generate correlated
rnum(1) = mean(1) + a11 * x1
@@ -147,4 +163,182 @@
!========================================================================
+
+!-------------------------------------------------------------------
+
+subroutine initialize_module
+
+if ( .not. module_initialized ) then
+ call register_module(source, revision, revdate)
+ module_initialized = .true.
+endif
+
+end subroutine initialize_module
+
+!-------------------------------------------------------------------
+
+! A random congruential random number generator of the Mersenne Twister type
+! See GNU Scientific Library code.
+
+subroutine init_ran(s, seed)
+ integer, intent(in) :: seed
+ type(random_seq_type), intent(out) :: s
+
+ integer :: i, sd
+
+if ( .not. module_initialized ) call initialize_module
+
+! Initialize the generator for use with
+! repeatable sequences
+
+sd = seed
+if (sd == 0) sd = 4357 ! do not allow seed to be 0, use default
+
+s%mt(0) = iand(sd, z'FFFFFFFF')
+
+! See Knuth's "Art of Computer Programming" Vol. 2, 3rd Ed. p.106
+! for multiplier.
+
+do i=1, N-1
+ s%mt(i) = 1812433253_i8 * (ieor(s%mt(i-1), ishft(s%mt(i-1), -30)) + i)
+
+ s%mt(i) = iand(s%mt(i), z'FFFFFFFF')
+end do
+
+s%mti = N
+s%gset = .false.
+
+end subroutine init_ran
+
+!-----------------------------------------------------------------
+
+! A random congruential random number generator from
+! the GNU Scientific Library (The Mersenne Twister MT19937 varient.)
+
+function ran_unif(s)
+ type(random_seq_type), intent(inout) :: s
+ real(r8) :: ran_unif
+
+integer :: kk
+integer(i8) :: k, y, m1
+
+if ( .not. module_initialized ) call initialize_module
+
+
+! 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
+ 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), 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-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-N)), ieor(ishft(y,-1), 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), 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))
+k = ieor(k, iand(ishft(k, 7), z'9d2c5680'))
+k = ieor(k, iand(ishft(k, 15), z'efc60000'))
+k = ieor(k, ishft(k, -18))
+
+! 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 wants a real between 0 and 1.0,
+! so divide here. i do not know if the expected
+! range is [0,1) (0,1] (0,1) or [0,1]. will search docs.
+
+ran_unif = real(k, r8) / 4294967296.0_r8
+
+end function ran_unif
+
+!------------------------------------------------------------------------
+
+function ran_gauss(s)
+
+! Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122
+
+! Returns a N(-1, 1) random number draw from a gaussian distribution
+
+type(random_seq_type), intent(inout) :: s
+real(r8) :: ran_gauss
+
+real(digits12) :: x, y, r2, t
+
+if ( .not. module_initialized ) call initialize_module
+
+10 continue
+
+! choose x,y in uniform square (-1,-1) to (+1,+1)
+x = -1.0_digits12 + 2.0_digits12 * ran_unif(s)
+y = -1.0_digits12 + 2.0_digits12 * ran_unif(s)
+
+! repeat if it is outside the unit circle or at origin
+r2 = x*x + y*y
+if (r2 > 1.0_digits12 .or. r2 == 0.0_digits12) goto 10
+
+if (s%gset) then
+ ran_gauss = s%lastg
+ s%gset = .false.
+else
+ t = sqrt(-2.0_digits12 * log(r2) / r2)
+ s%lastg = x * t
+ s%gset = .true.
+ ran_gauss = y * t
+endif
+
+end function ran_gauss
+
+!------------------------------------------------------------------------
+
+
end module random_seq_mod
Modified: DART/branches/development/random_seq/random_seq_mod.html
===================================================================
--- DART/branches/development/random_seq/random_seq_mod.html 2012-03-16 23:11:46 UTC (rev 5608)
+++ DART/branches/development/random_seq/random_seq_mod.html 2012-03-19 14:50:48 UTC (rev 5609)
@@ -50,7 +50,6 @@
<PRE>
types_mod
utilities_mod
-random_nr_mod
</PRE>
<!--==================================================================-->
@@ -66,7 +65,7 @@
<H2>PUBLIC INTERFACES</H2>
<TABLE>
-<TR><TD><em class=call>use assim_model_mod, only : </em></TD>
+<TR><TD><em class=call>use random_seq_mod, only : </em></TD>
<TD><A HREF="#random_seq_type">random_seq_type</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#init_random_seq">init_random_seq</A></TD></TR>
<TR><TD> </TD><TD><A HREF="#random_gaussian">random_gaussian</A></TD></TR>
@@ -95,6 +94,12 @@
<pre>
<em class=call>type random_seq_type</em>
private
+
+ integer :: mti
+ integer(i8) :: mt(624)
+ real(r8) :: lastg
+ logical :: gset
+
end type random_seq_type
</pre>
</div>
@@ -103,9 +108,9 @@
<!-- Description -->
<P>
-This type is used to uniquely identify a sequence. In this implementation,
-it is passed through directly from the random_nr_mod and details can be
-seen in that module.
+This type is used to uniquely identify a sequence. Keeps the state history
+of the linear congruential number generator. In this implementation it is
+based on the Mersenne Twister from the GNU Scientific Library.
</P>
</div>
@@ -132,7 +137,6 @@
of independent, reproducible random sequences can be generated by
having multiple instances of a random_seq_type. A specified integer
seed, optional, can produce a specific 'random' sequence.
-Seed must be a negative number.
</P>
<TABLE width=100% border=0 summary="" cellpadding=3>
@@ -314,7 +318,14 @@
<A NAME="References"></A>
<HR>
<H2>REFERENCES</H2>
+<ol>
+<li> Knuth, Vol 2. </li>
+<li> <a href="http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html">
+GNU Scientific Library Reference Manual</a> </li>
+</li>
+</ol>
+
<!--==================================================================-->
<!-- Describe all the error conditions and codes. -->
<!--==================================================================-->
@@ -356,10 +367,103 @@
<A NAME="PrivateComponents"></A>
<HR>
<H2>PRIVATE COMPONENTS</H2>
+
+<TABLE>
+<TR><TD> </TD><TD><A HREF="#init_ran">init_ran</A></TD></TR>
+<TR><TD> </TD><TD><A HREF="#ran_unif">ran_unif</A></TD></TR>
+<TR><TD> </TD><TD><A HREF="#ran_gauss">ran_gauss</A></TD></TR>
+</TABLE>
+
+
+<!--===================== DESCRIPTION OF A ROUTINE =====================-->
+
+<A NAME="init_ran"></A>
+<br>
+<div class=routine>
+<em class=call> call init_ran(s,temp) </em>
+<pre>
+type(random_seq_type), intent(out) :: <em class=code>s</em>
+integer, intent(in) :: <em class=code>temp</em>
+</pre>
+</div>
+
+<div class=indent1>
+<!-- Description -->
+
<P>
-N/A
+Initializes a random sequence with an integer. Any sequence
+initialized with the same integer will produce the same sequence
+of pseudo-random numbers.
</P>
+<TABLE width=100% border=0 summary="" cellpadding=3>
+<TR><TD valign=top><em class=code>s </em></TD>
+ <TD>A random sequence to be initialized</TD></TR>
+<TR><TD valign=top><em class=code>temp </em></TD>
+ <TD>An integer seed to start the sequence.</TD></TR>
+</TABLE>
+
+</div>
+<br>
+
+<!--===================== DESCRIPTION OF A ROUTINE =====================-->
+
+<A NAME="ran_unif"></A>
+<br>
+<div class=routine>
+<em class=call> var = ran_unif(s) </em>
+<pre>
+real(r8) :: <em class=code>ran_unif</em>
+type(random_seq_type), intent(inout) :: <em class=code>s</em>
+</pre>
+</div>
+
+<div class=indent1>
+<!-- Description -->
+
+<P>
+Generate the next uniform [0, 1] random number in the sequence.
+</P>
+
+<TABLE width=100% border=0 summary="" cellpadding=3>
+<TR><TD valign=top><em class=code>ran1 </em></TD>
+ <TD>Next uniformly distributed [0, 1] number in sequence.</TD></TR>
+<TR><TD valign=top><em class=code>s </em></TD>
+ <TD>A random sequence.</TD></TR>
+</TABLE>
+
+</div>
+<br>
+
+<!--===================== DESCRIPTION OF A ROUTINE =====================-->
+
+<A NAME="ran_gauss"></A>
+<br>
+<div class=routine>
+<em class=call> var = ran_gauss(s) </em>
+<pre>
+real(r8) :: <em class=code>ran_gauss</em>
+type(random_seq_type), intent(inout) :: <em class=code>s</em>
+</pre>
+</div>
+
+<div class=indent1>
+<!-- Description -->
+
+<P>
+Generates a random draw from a standard gaussian.
+</P>
+
+<TABLE width=100% border=0 summary="" cellpadding=3>
+<TR><TD valign=top><em class=code>ran_gauss </em></TD>
+ <TD>A random draw from a standard gaussian.</TD></TR>
+<TR><TD valign=top><em class=code>s </em></TD>
+ <TD>A random sequence.</TD></TR>
+</TABLE>
+
+</div>
+<br>
+
<!--==================================================================-->
<!-- Legalese & Metadata -->
<!--==================================================================-->
Added: DART/branches/development/random_seq/test/input.nml
===================================================================
--- DART/branches/development/random_seq/test/input.nml (rev 0)
+++ DART/branches/development/random_seq/test/input.nml 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,20 @@
+&utilities_nml
+ TERMLEVEL = 1,
+ module_details = .false.
+ logfilename = 'dart_log.out'
+/
+
+&test_diff_nml
+ mean = 0.0d0,
+ sd1 = 1.0d0,
+ sd2 = 1.0d0,
+ iseed = 123
+/
+
+&test_reseed_nml
+ nsamples = 1000
+ calendar_name = 'NO_CALENDAR'
+ start_day = 0
+ start_sec = 0
+/
+
Property changes on: DART/branches/development/random_seq/test/input.nml
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:eol-style
+ native
Added: DART/branches/development/random_seq/test/mkmf_test_corr
===================================================================
--- DART/branches/development/random_seq/test/mkmf_test_corr (rev 0)
+++ DART/branches/development/random_seq/test/mkmf_test_corr 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+../../mkmf/mkmf -p test_corr -t ../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../" path_names_test_corr
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/branches/development/random_seq/test/mkmf_test_corr
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Copied: DART/branches/development/random_seq/test/mkmf_test_diff (from rev 5592, DART/branches/development/random_seq/mkmf_test_diff)
===================================================================
--- DART/branches/development/random_seq/test/mkmf_test_diff (rev 0)
+++ DART/branches/development/random_seq/test/mkmf_test_diff 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+../../mkmf/mkmf -p test_diff -t ../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../" path_names_test_diff
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Added: DART/branches/development/random_seq/test/mkmf_test_random_gsl
===================================================================
--- DART/branches/development/random_seq/test/mkmf_test_random_gsl (rev 0)
+++ DART/branches/development/random_seq/test/mkmf_test_random_gsl 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+../../mkmf/mkmf -p test_random_gsl -t ../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../" path_names_test_random_gsl
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/branches/development/random_seq/test/mkmf_test_random_gsl
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/random_seq/test/mkmf_test_reseed
===================================================================
--- DART/branches/development/random_seq/test/mkmf_test_reseed (rev 0)
+++ DART/branches/development/random_seq/test/mkmf_test_reseed 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,18 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+../../mkmf/mkmf -p test_reseed -t ../../mkmf/mkmf.template -c"-Duse_netCDF" \
+ -a "../../" path_names_test_reseed
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/branches/development/random_seq/test/mkmf_test_reseed
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/random_seq/test/path_names_test_corr
===================================================================
--- DART/branches/development/random_seq/test/path_names_test_corr (rev 0)
+++ DART/branches/development/random_seq/test/path_names_test_corr 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,6 @@
+random_seq/test_corr.f90
+random_seq/random_seq_mod.f90
+common/types_mod.f90
+utilities/utilities_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+time_manager/time_manager_mod.f90
Copied: DART/branches/development/random_seq/test/path_names_test_diff (from rev 5592, DART/branches/development/random_seq/path_names_test_diff)
===================================================================
--- DART/branches/development/random_seq/test/path_names_test_diff (rev 0)
+++ DART/branches/development/random_seq/test/path_names_test_diff 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,6 @@
+random_seq/test_diff.f90
+random_seq/random_seq_mod.f90
+common/types_mod.f90
+utilities/utilities_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+time_manager/time_manager_mod.f90
Added: DART/branches/development/random_seq/test/path_names_test_random_gsl
===================================================================
--- DART/branches/development/random_seq/test/path_names_test_random_gsl (rev 0)
+++ DART/branches/development/random_seq/test/path_names_test_random_gsl 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,6 @@
+random_seq/test_random_gsl.f90
+random_seq/random_seq_mod.f90
+common/types_mod.f90
+utilities/utilities_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+time_manager/time_manager_mod.f90
Added: DART/branches/development/random_seq/test/path_names_test_reseed
===================================================================
--- DART/branches/development/random_seq/test/path_names_test_reseed (rev 0)
+++ DART/branches/development/random_seq/test/path_names_test_reseed 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,6 @@
+random_seq/test_reseed.f90
+random_seq/random_seq_mod.f90
+common/types_mod.f90
+utilities/utilities_mod.f90
+mpi_utilities/null_mpi_utilities_mod.f90
+time_manager/time_manager_mod.f90
Added: DART/branches/development/random_seq/test_random_gsl.f90
===================================================================
--- DART/branches/development/random_seq/test_random_gsl.f90 (rev 0)
+++ DART/branches/development/random_seq/test_random_gsl.f90 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,86 @@
+! DART software - Copyright 2004 - 2011 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_random_gsl
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+use types_mod, only : r4, r8, digits12
+use utilities_mod, only : register_module, error_handler, E_ERR, &
+ initialize_utilities, timestamp
+use random_mod, only : random_seq_type, init_ran, ran_unif, ran_gauss
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+type (random_seq_type) :: r
+integer :: i, n
+
+double precision :: dpr1, dpdist, dpmean_dist
+real(r8) :: r8r1, r8dist, r8mean_dist
+real(r4) :: r4r1, r4dist, r4mean_dist
+real(digits12) :: d12r1, d12dist, d12mean_dist
+
+call initialize_utilities('test_random')
+call register_module(source,revision,revdate)
+
+write(*, *) 'double precision is ', kind(dpr1)
+write(*, *) 'digits12 is defined to be ', digits12
+write(*, *) 'r8 is defined to be ', r8
+write(*, *) 'r4 is defined to be ', r4
+
+n = 10000000
+
+
+call init_ran(r, -5)
+d12mean_dist = 0.0
+do i = 1, n
+ d12r1 = ran_gauss(r)
+ d12dist = abs(d12r1)
+ d12mean_dist = d12mean_dist + d12dist
+end do
+write(*, *) 'digits12 sd is ', d12mean_dist / n
+
+
+call init_ran(r, -5)
+dpmean_dist = 0.0
+do i = 1, n
+ dpr1 = ran_gauss(r)
+ dpdist = dabs(dpr1)
+ dpmean_dist = dpmean_dist + dpdist
+end do
+write(*, *) 'double precision sd is ', dpmean_dist / n
+
+
+call init_ran(r, -5)
+r8mean_dist = 0.0_r8
+do i = 1, n
+ r8r1 = ran_gauss(r)
+ r8dist = abs(r8r1)
+ r8mean_dist = r8mean_dist + r8dist
+end do
+write(*, *) 'r8 sd is ', r8mean_dist / n
+
+
+call init_ran(r, -5)
+r4mean_dist = 0.0_r4
+do i = 1, n
+ r4r1 = ran_gauss(r)
+ r4dist = abs(r4r1)
+ r4mean_dist = r4mean_dist + r4dist
+end do
+write(*, *) 'r4 sd is ', r4mean_dist / n
+
+call timestamp(source,revision,revdate,'end')
+
+end program test_random_gsl
Property changes on: DART/branches/development/random_seq/test_random_gsl.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/random_seq/test_reseed.f90
===================================================================
--- DART/branches/development/random_seq/test_reseed.f90 (rev 0)
+++ DART/branches/development/random_seq/test_reseed.f90 2012-03-19 14:50:48 UTC (rev 5609)
@@ -0,0 +1,340 @@
+! DART software - Copyright 2004 - 2011 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_reseed
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+use types_mod, only : r4, r8, digits12, i8
+use utilities_mod, only : register_module, error_handler, E_ERR, &
+ initialize_utilities, finalize_utilities, &
+ logfileunit, nmlfileunit, timestamp, &
+ find_namelist_in_file, check_namelist_read, &
+ open_file, close_file, do_nml_file, do_nml_term
+use time_manager_mod, only : time_type, operator(+), set_time, get_time, &
+ set_calendar_type, print_time, print_date
+use random_seq_mod, only : random_seq_type, init_random_seq, &
+ random_uniform, random_gaussian
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+type (time_type) :: state_time
+
+integer :: io
+integer :: unitnum = 66, iunit
+character(len=128) :: fname
+integer :: reseed, hours
+
+character(len=64) :: calendar_name = 'GREGORIAN'
+integer :: nsamples = 1000000
+integer :: start_day = 148866
+integer :: start_sec = 0
+
+namelist /test_reseed_nml/ calendar_name, nsamples, &
+start_day, start_sec
+
+! --------------------------------------------
+
+call initialize_utilities('test_reseed')
+call register_module(source,revision,revdate)
+
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "test_reseed_nml", iunit)
+read(iunit, nml = test_reseed_nml, iostat = io)
+call check_namelist_read(iunit, io, "test_reseed_nml")
+
+! Record the namelist values used for the run ...
+if (do_nml_file()) write(nmlfileunit, nml=test_reseed_nml)
+if (do_nml_term()) write( * , nml=test_reseed_nml)
+
+! test no calendar like low order models use, also test
+! the gregorian cal with a date close to current.
+
+call set_calendar_type(calendar_name)
+state_time = set_time(start_sec, start_day)
+
+print *, ' '
+call print_time(state_time, 'setting start time')
+if (calendar_name == 'GREGORIAN') call print_date(state_time, 'is start date')
+
+! -----
+
+fname = 'uniform_baseline'
+open(unit=unitnum, file=fname, action='write')
+call test1(nsamples, .true., state_time, 66)
+close(unitnum)
+
+fname = 'gaussian_baseline'
+open(unit=unitnum, file=fname, action='write')
+call test1(nsamples, .false., state_time, 66)
+close(unitnum)
+
+! -----
+
+hours = 1
+reseed = 1
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 1
+reseed = 10
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 6
+reseed = 10
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 12
+reseed = 10
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 6
+reseed = 100
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 12
+reseed = 100
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 24
+reseed = 100
+
+call dohourstest(reseed, hours)
+
+! -----
+
+hours = 24
+reseed = 1000
+
+call dohourstest(reseed, hours)
+
+
+! -----
+
+call finalize_utilities()
+
+contains
+
+!-------------------------------------------------------------------------
+
+subroutine test1(nreps, unif, start_time, unit)
+! output uniform distribution. seeds gen with a single time
+ integer, intent(in) :: nreps, unit
+ logical, intent(in) :: unif
+ type(time_type), intent(in) :: start_time
+
+integer :: i, seed
+type (random_seq_type) :: seq1
+
+seed = generate_seed(start_time)
+call init_random_seq(seq1, seed)
+
+do i=1, nreps
+ if (unif) then
+ write(unit, *) random_uniform(seq1)
+ else
+ write(unit, *) random_gaussian(seq1, 0.0_r8, 1.0_r8)
+ endif
+end do
+
+end subroutine test1
+
+!-------------------------------------------------------------------------
+
+subroutine test2(nreps, per, start_time, delta_time, unif, unit)
+! output uniform distribution. generates nrep numbers, reseeding
+! the generator each 'per' times.
+ integer, intent(in) :: nreps, per, unit
+ type(time_type), intent(in) :: start_time, delta_time
+ logical, intent(in) :: unif
+
+integer :: i, j, nextseed
+type(time_type) :: t
+type(random_seq_type) :: seq1
+
+t = start_time
+nextseed = generate_seed(t)
+call init_random_seq(seq1, nextseed)
+
+j = 1
+do i=1, nreps
+ if (unif) then
+ write(unit, *) random_uniform(seq1)
+ else
+ write(unit, *) random_gaussian(seq1, 0.0_r8, 1.0_r8)
+ endif
+ if (j == per) then
+ t = t + delta_time
+ nextseed = generate_seed(t)
+ call init_random_seq(seq1, nextseed)
+ j = 1
+ else
+ j = j + 1
+ endif
+end do
+
+end subroutine test2
+
+!-------------------------------------------------------------------------
+
+subroutine test3
+
+real(r8), allocatable :: history(:)
+real(r8) :: next_val
+integer :: i, j, nextseed
+type(time_type) :: t, base_time, state_time, delta_time, delta_time2
+type(random_seq_type) :: seq1
+
+
+base_time = set_time(0, 148000)
+state_time = base_time
+delta_time = set_time(1)
+delta_time2 = set_time(0, 1)
+
+print *, ' '
+call print_time(state_time, 'setting start time')
+call print_time(delta_time, 'delta time')
+call print_time(delta_time2, 'delta time2')
+print *, ' '
+
+! generate the first number from 100 consecutive seeds
+allocate(history(10000))
+do i=1, 10000
+ nextseed = generate_seed(state_time)
+
+ call init_random_seq(seq1, nextseed)
+ history(i) = random_uniform(seq1)
+
+ state_time = state_time + delta_time
+enddo
+
+! now generate 100 values from each seed and see if
+! they overlap in any of the previously generated vals.
+state_time = base_time + delta_time2
+do i=1, 1000
+ nextseed = generate_seed(state_time)
+
+ call init_random_seq(seq1, nextseed)
+ next_val = random_uniform(seq1)
+
+ do j=1, 10000-1
+ if (history(j) == next_val) then
+ print *, 'found match, ', i, j, history(j), next_val
+ next_val = random_uniform(seq1)
+ print *, 'next val from ran ', next_val, ' next in seq is ', history(j+1)
+ endif
+ enddo
+
+ state_time = state_time + delta_time2
+
+enddo
+
+deallocate(history)
+
+end subroutine test3
+
+!-------------------------------------------------------------------------
+
+subroutine dohourstest(reseed, hours)
+ integer, intent(in) :: reseed, hours
+
+type(time_type) :: delta_time
+
+write(fname, '(A,I4.4,A,I2.2,A)') 'uniform_', reseed, 'dt_', hours, 'h'
+open(unit=unitnum, file=fname, action='write')
+delta_time = set_time(hours*3600, 0)
+call test2(nsamples, reseed, state_time, delta_time, .true., unitnum)
+close(unitnum)
+
+write(fname, '(A,I4.4,A,I2.2,A)') 'gaussian_', reseed, 'dt_', hours, 'h'
+open(unit=unitnum, file=fname, action='write')
+delta_time = set_time(hours*3600, 0)
+call test2(nsamples, reseed, state_time, delta_time, .false., unitnum)
+close(unitnum)
+
+end subroutine dohourstest
+
+!-------------------------------------------------------------------------
+
+function generate_seed(state_time)
+! use the state time to set the seed for the (repeatable) random sequence
+
+type(time_type), intent(in) :: state_time
+integer :: generate_seed
+
+integer :: days,seconds,bigint,bigtwo
+integer(i8) :: bigtime,bigone
+integer(i8), parameter :: secs_day = 86400_i8
+
+call get_time(state_time, seconds, days)
+
+! option 1: add days & secs, old gen needs negative seed,
+! new one doesn't care.
+!generate_seed = days + seconds
+!return
+
+! option 2: compute total seconds in an i8 then use the lower
+! 32 bits (which is all the init routine will take to stay
+! portable across different platforms).
+
+! if generate_seed was i8, or if this routine was combined with
+! the real init routine, we could avoid the bitwise and below.
+! it's going to get repeated in the set seed routine.
+generate_seed = iand((secs_day * days) + seconds, z'FFFFFFFF')
+return
+
+! option 3: works with old generator but maybe overkill?
+! for the old generator, it has to be < 0 and apparently > something
+! like 50,000 more than the largest negative int.
+
+!bigtime = int(days,i8)*100000_i8 + int(seconds,i8)
+!bigint = huge(seconds) ! biggest 32bit integer
+!bigone = int(bigint,i8) ! coerce to 64bit integer
+!bigtwo = int(mod(bigtime,bigone)) ! modulo arith on 64bit integers, then 32bit
+!generate_seed = -1 * int(bigtwo) ! coerce back to 32bit integer
+
+! arb choices if the seed is going to cause an error
+!if (generate_seed >= 0) generate_seed = -1234
+!if (generate_seed < -2147432706) generate_seed = -4376
+!return
+
+! option 4: positive number instead of negative
+!generate_seed = bigtwo
+!return
+
+! not reached anymore:
+!write(*,*)'TJH DEBUG generate_seed ',bigtime,bigint,generate_seed
+
+end function generate_seed
+
+!-------------------------------------------------------------------------
+
+end program test_reseed
Property changes on: DART/branches/development/random_seq/test_reseed.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
More information about the Dart-dev
mailing list