(git:374b731)
Loading...
Searching...
No Matches
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
8PROGRAM 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
132CONTAINS
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
245END 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::advance
subroutine advance(self, e, c)
Advance the state by n steps, i.e. jump n steps forward, if n > 0, or backward if n < 0.
Definition parallel_rng_types.F:228
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
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
message_passing::mp_comm_type
Definition message_passing.F:189
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157