(git:e7e05ae)
parallel_rng_types_unittest.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 PROGRAM parallel_rng_types_test
9  USE message_passing, ONLY: mp_world_finalize, &
10  mp_world_init, &
11  mp_comm_type
12  USE kinds, ONLY: dp
13  USE machine, ONLY: m_walltime, &
14  default_output_unit
15  USE parallel_rng_types, ONLY: gaussian, &
16  uniform, &
17  check_rng, &
18  rng_stream_type, &
19  rng_stream_type_from_record, &
20  rng_name_length, &
21  rng_record_length
22 
23  IMPLICIT NONE
24 
25  INTEGER :: i, nsamples, nargs, stat
26  LOGICAL :: ionode
27  REAL(kind=dp) :: t, tend, tmax, tmin, tstart, tsum, tsum2
28  TYPE(mp_comm_type) :: mpi_comm
29  TYPE(rng_stream_type) :: rng_stream
30  CHARACTER(len=32) :: arg
31 
32  nsamples = 1000
33  nargs = command_argument_count()
34 
35  IF (nargs .GT. 1) &
36  error stop "Usage: parallel_rng_types_TEST [<int:nsamples>]"
37 
38  IF (nargs == 1) THEN
39  CALL get_command_argument(1, arg)
40  READ (arg, *, iostat=stat) nsamples
41  IF (stat /= 0) &
42  error stop "Usage: parallel_rng_types_TEST [<int:nsamples>]"
43  END IF
44 
45  CALL mp_world_init(mpi_comm)
46  ionode = mpi_comm%is_source()
47 
48  CALL check_rng(default_output_unit, ionode)
49 
50  ! Check performance
51 
52  IF (ionode) THEN
53  WRITE (unit=default_output_unit, fmt="(/,/,T2,A,I10,A)") &
54  "Check distributions using", nsamples, " random numbers:"
55  END IF
56 
57  ! Test uniform distribution [0,1]
58 
59  rng_stream = rng_stream_type(name="Test uniform distribution [0,1]", &
60  distribution_type=uniform, &
61  extended_precision=.true.)
62 
63  IF (ionode) &
64  CALL rng_stream%write(default_output_unit)
65 
66  tmax = -huge(0.0_dp)
67  tmin = +huge(0.0_dp)
68  tsum = 0.0_dp
69  tsum2 = 0.0_dp
70 
71  tstart = m_walltime()
72  DO i = 1, nsamples
73  t = rng_stream%next()
74  tsum = tsum + t
75  tsum2 = tsum2 + t*t
76  IF (t > tmax) tmax = t
77  IF (t < tmin) tmin = t
78  END DO
79  tend = m_walltime()
80 
81  IF (ionode) THEN
82  CALL rng_stream%write(default_output_unit, write_all=.true.)
83  WRITE (unit=default_output_unit, fmt="(/,(T4,A,F12.6))") &
84  "Minimum: ", tmin, &
85  "Maximum: ", tmax, &
86  "Average: ", tsum/real(nsamples, kind=dp), &
87  "Variance:", tsum2/real(nsamples, kind=dp), &
88  "Time [s]:", tend - tstart
89  END IF
90 
91  ! Test normal Gaussian distribution
92 
93  rng_stream = rng_stream_type(name="Test normal Gaussian distribution", &
94  distribution_type=gaussian, &
95  extended_precision=.true.)
96 
97  IF (ionode) &
98  CALL rng_stream%write(default_output_unit)
99 
100  tmax = -huge(0.0_dp)
101  tmin = +huge(0.0_dp)
102  tsum = 0.0_dp
103  tsum2 = 0.0_dp
104 
105  tstart = m_walltime()
106  DO i = 1, nsamples
107  t = rng_stream%next()
108  tsum = tsum + t
109  tsum2 = tsum2 + t*t
110  IF (t > tmax) tmax = t
111  IF (t < tmin) tmin = t
112  END DO
113  tend = m_walltime()
114 
115  IF (ionode) THEN
116  CALL rng_stream%write(default_output_unit)
117  WRITE (unit=default_output_unit, fmt="(/,(T4,A,F12.6))") &
118  "Minimum: ", tmin, &
119  "Maximum: ", tmax, &
120  "Average: ", tsum/real(nsamples, kind=dp), &
121  "Variance:", tsum2/real(nsamples, kind=dp), &
122  "Time [s]:", tend - tstart
123  END IF
124 
125  IF (ionode) THEN
126  CALL dump_reload_check()
127  CALL shuffle_check()
128  END IF
129 
130  CALL mp_world_finalize()
131 
132 CONTAINS
133 ! **************************************************************************************************
134 !> \brief ...
135 ! **************************************************************************************************
136  SUBROUTINE dump_reload_check()
137  TYPE(rng_stream_type) :: rng_stream
138  CHARACTER(len=rng_record_length) :: rng_record
139  REAL(KIND=dp), DIMENSION(3, 2) :: ig, ig_orig, cg, cg_orig, bg, bg_orig
140  CHARACTER(len=rng_name_length) :: name, name_orig
141  CHARACTER(len=*), PARAMETER :: serialized_string = &
142  "qtb_rng_gaussian 1 F T F 0.0000000000000000E+00&
143  & 12.0 12.0 12.0&
144  & 12.0 12.0 12.0&
145  & 12.0 12.0 12.0&
146  & 12.0 12.0 12.0&
147  & 12.0 12.0 12.0&
148  & 12.0 12.0 12.0"
149 
150  WRITE (unit=default_output_unit, fmt="(/,/,T2,A)") &
151  "Checking dump and load round trip:"
152 
153  rng_stream = rng_stream_type(name="Roundtrip for normal Gaussian distrib", &
154  distribution_type=gaussian, &
155  extended_precision=.true.)
156 
157  CALL rng_stream%advance(7, 42)
158  CALL rng_stream%get(ig=ig_orig, cg=cg_orig, bg=bg_orig, name=name_orig)
159  CALL rng_stream%dump(rng_record)
160 
161  rng_stream = rng_stream_type_from_record(rng_record)
162  CALL rng_stream%get(ig=ig, cg=cg, bg=bg, name=name)
163 
164  IF (any(ig /= ig_orig) .OR. any(cg /= cg_orig) .OR. any(bg /= bg_orig) &
165  .OR. (name /= name_orig)) &
166  error stop "Stream dump and load roundtrip failed"
167 
168  WRITE (unit=default_output_unit, fmt="(T4,A)") &
169  "Roundtrip successful"
170 
171  WRITE (unit=default_output_unit, fmt="(/,/,T2,A)") &
172  "Checking dumped format:"
173 
174  ig(:, :) = 12.0_dp
175  rng_stream = rng_stream_type(name="qtb_rng_gaussian", &
176  distribution_type=gaussian, &
177  extended_precision=.true., &
178  seed=ig)
179 
180  CALL rng_stream%dump(rng_record)
181 
182  WRITE (unit=default_output_unit, fmt="(T4,A10,A433)") &
183  "EXPECTED:", serialized_string
184 
185  WRITE (unit=default_output_unit, fmt="(T4,A10,A433)") &
186  "GENERATED:", rng_record
187 
188  IF (rng_record /= serialized_string) &
189  error stop "Serialized record does not match the expected output"
190 
191  WRITE (unit=default_output_unit, fmt="(T4,A)") &
192  "Serialized record matches the expected output"
193 
194  END SUBROUTINE
195 
196 ! **************************************************************************************************
197 !> \brief ...
198 ! **************************************************************************************************
199  SUBROUTINE shuffle_check()
200  TYPE(rng_stream_type) :: rng_stream
201 
202  INTEGER, PARAMETER :: sz = 20
203  INTEGER, DIMENSION(1:sz) :: arr, arr2, orig
204  LOGICAL, DIMENSION(1:sz) :: mask
205  INTEGER :: idx
206  REAL(KIND=dp), DIMENSION(3, 2), PARAMETER :: ig = 12.0_dp
207 
208  WRITE (unit=default_output_unit, fmt="(/,/,T2,A)", advance="no") &
209  "Checking shuffle()"
210 
211  rng_stream = rng_stream_type(name="shuffle() check", seed=ig)
212  orig = [(idx, idx=1, sz)]
213 
214  arr = orig
215  CALL rng_stream%shuffle(arr)
216 
217  IF (all(arr == orig)) &
218  error stop "shuffle failed: array was left untouched"
219  WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
220 
221  IF (any(arr /= orig(arr))) &
222  error stop "shuffle failed: the shuffled original is not the shuffled original"
223  WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
224 
225  ! sort and compare to orig
226  mask = .true.
227  DO idx = 1, size(orig)
228  IF (minval(arr, mask) /= orig(idx)) &
229  error stop "shuffle failed: there is at least one unknown index"
230  mask(minloc(arr, mask)) = .false.
231  END DO
232  WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
233 
234  arr2 = orig
235  CALL rng_stream%reset()
236  CALL rng_stream%shuffle(arr2)
237 
238  IF (any(arr2 /= arr)) &
239  error stop "shuffle failed: array was shuffled differently with same rng state"
240  WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
241 
242  WRITE (unit=default_output_unit, fmt="(T4,A)") &
243  " successful"
244  END SUBROUTINE
245 END PROGRAM parallel_rng_types_test
246 ! vim: set ts=3 sw=3 tw=132 :
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition: machine.F:45
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
message_passing::mp_world_init
subroutine, public mp_world_init(mp_comm)
initializes the system default communicator
Definition: message_passing.F:1274
message_passing::mp_world_finalize
subroutine, public mp_world_finalize()
finalizes the system default communicator
Definition: message_passing.F:1386
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition: parallel_rng_types.F:51
parallel_rng_types::rng_stream_type_from_record
type(rng_stream_type) function, public rng_stream_type_from_record(rng_record)
Create a RNG stream from a record given as an internal file (string).
Definition: parallel_rng_types.F:564
parallel_rng_types::rng_name_length
integer, parameter, public rng_name_length
Definition: parallel_rng_types.F:66
parallel_rng_types::rng_record_length
integer, parameter, public rng_record_length
Definition: parallel_rng_types.F:65
parallel_rng_types::uniform
integer, parameter, public uniform
Definition: parallel_rng_types.F:73
parallel_rng_types::check_rng
subroutine, public check_rng(output_unit, ionode)
...
Definition: parallel_rng_types.F:279
parallel_rng_types::gaussian
integer, parameter, public gaussian
Definition: parallel_rng_types.F:73
shuffle_check
subroutine shuffle_check()
...
Definition: parallel_rng_types_unittest.F:200
dump_reload_check
subroutine dump_reload_check()
...
Definition: parallel_rng_types_unittest.F:137
parallel_rng_types_test
program parallel_rng_types_test
Definition: parallel_rng_types_unittest.F:8