(git:d2a9ebd)
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-2026 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 > 1) then
36 error stop "Usage: parallel_rng_types_TEST [<int:nsamples>]"
37 end if
38
39 IF (nargs == 1) THEN
40 CALL get_command_argument(1, arg)
41 READ (arg, *, iostat=stat) nsamples
42 IF (stat /= 0) then
43 error stop "Usage: parallel_rng_types_TEST [<int:nsamples>]"
44 end if
45 END IF
46
47 CALL mp_world_init(mpi_comm)
48 ionode = mpi_comm%is_source()
49
50 CALL check_rng(default_output_unit, ionode)
51
52 ! Check performance
53
54 IF (ionode) THEN
55 WRITE (unit=default_output_unit, fmt="(/,/,T2,A,I10,A)") &
56 "Check distributions using", nsamples, " random numbers:"
57 END IF
58
59 ! Test uniform distribution [0,1]
60
61 rng_stream = rng_stream_type(name="Test uniform distribution [0,1]", &
62 distribution_type=uniform, &
63 extended_precision=.true.)
64
65 IF (ionode) then
66 CALL rng_stream%write(default_output_unit)
67 end if
68
69 tmax = -huge(0.0_dp)
70 tmin = +huge(0.0_dp)
71 tsum = 0.0_dp
72 tsum2 = 0.0_dp
73
74 tstart = m_walltime()
75 DO i = 1, nsamples
76 t = rng_stream%next()
77 tsum = tsum + t
78 tsum2 = tsum2 + t*t
79 IF (t > tmax) tmax = t
80 IF (t < tmin) tmin = t
81 END DO
82 tend = m_walltime()
83
84 IF (ionode) THEN
85 CALL rng_stream%write(default_output_unit, write_all=.true.)
86 WRITE (unit=default_output_unit, fmt="(/,(T4,A,F12.6))") &
87 "Minimum: ", tmin, &
88 "Maximum: ", tmax, &
89 "Average: ", tsum/real(nsamples, kind=dp), &
90 "Variance:", tsum2/real(nsamples, kind=dp), &
91 "Time [s]:", tend - tstart
92 END IF
93
94 ! Test normal Gaussian distribution
95
96 rng_stream = rng_stream_type(name="Test normal Gaussian distribution", &
97 distribution_type=gaussian, &
98 extended_precision=.true.)
99
100 IF (ionode) then
101 CALL rng_stream%write(default_output_unit)
102 end if
103
104 tmax = -huge(0.0_dp)
105 tmin = +huge(0.0_dp)
106 tsum = 0.0_dp
107 tsum2 = 0.0_dp
108
109 tstart = m_walltime()
110 DO i = 1, nsamples
111 t = rng_stream%next()
112 tsum = tsum + t
113 tsum2 = tsum2 + t*t
114 IF (t > tmax) tmax = t
115 IF (t < tmin) tmin = t
116 END DO
117 tend = m_walltime()
118
119 IF (ionode) THEN
120 CALL rng_stream%write(default_output_unit)
121 WRITE (unit=default_output_unit, fmt="(/,(T4,A,F12.6))") &
122 "Minimum: ", tmin, &
123 "Maximum: ", tmax, &
124 "Average: ", tsum/real(nsamples, kind=dp), &
125 "Variance:", tsum2/real(nsamples, kind=dp), &
126 "Time [s]:", tend - tstart
127 END IF
128
129 IF (ionode) THEN
130 CALL dump_reload_check()
131 CALL shuffle_check()
132 END IF
133
134 CALL mp_world_finalize()
135
136CONTAINS
137! **************************************************************************************************
138!> \brief ...
139! **************************************************************************************************
140 SUBROUTINE dump_reload_check()
141 TYPE(rng_stream_type) :: rng_stream
142 CHARACTER(len=rng_record_length) :: rng_record
143 REAL(KIND=dp), DIMENSION(3, 2) :: ig, ig_orig, cg, cg_orig, bg, bg_orig
144 CHARACTER(len=rng_name_length) :: name, name_orig
145 CHARACTER(len=*), PARAMETER :: serialized_string = &
146 "qtb_rng_gaussian 1 F T F 0.0000000000000000E+00&
147 & 12.0 12.0 12.0&
148 & 12.0 12.0 12.0&
149 & 12.0 12.0 12.0&
150 & 12.0 12.0 12.0&
151 & 12.0 12.0 12.0&
152 & 12.0 12.0 12.0"
153
154 WRITE (unit=default_output_unit, fmt="(/,/,T2,A)") &
155 "Checking dump and load round trip:"
156
157 rng_stream = rng_stream_type(name="Roundtrip for normal Gaussian distrib", &
158 distribution_type=gaussian, &
159 extended_precision=.true.)
160
161 CALL rng_stream%advance(7, 42)
162 CALL rng_stream%get(ig=ig_orig, cg=cg_orig, bg=bg_orig, name=name_orig)
163 CALL rng_stream%dump(rng_record)
164
165 rng_stream = rng_stream_type_from_record(rng_record)
166 CALL rng_stream%get(ig=ig, cg=cg, bg=bg, name=name)
167
168 IF (any(ig /= ig_orig) .OR. any(cg /= cg_orig) .OR. any(bg /= bg_orig) &
169 .OR. (name /= name_orig)) then
170 error stop "Stream dump and load roundtrip failed"
171 end if
172
173 WRITE (unit=default_output_unit, fmt="(T4,A)") &
174 "Roundtrip successful"
175
176 WRITE (unit=default_output_unit, fmt="(/,/,T2,A)") &
177 "Checking dumped format:"
178
179 ig(:, :) = 12.0_dp
180 rng_stream = rng_stream_type(name="qtb_rng_gaussian", &
181 distribution_type=gaussian, &
182 extended_precision=.true., &
183 seed=ig)
184
185 CALL rng_stream%dump(rng_record)
186
187 WRITE (unit=default_output_unit, fmt="(T4,A10,A433)") &
188 "EXPECTED:", serialized_string
189
190 WRITE (unit=default_output_unit, fmt="(T4,A10,A433)") &
191 "GENERATED:", rng_record
192
193 IF (rng_record /= serialized_string) then
194 error stop "Serialized record does not match the expected output"
195 end if
196
197 WRITE (unit=default_output_unit, fmt="(T4,A)") &
198 "Serialized record matches the expected output"
199
200 END SUBROUTINE dump_reload_check
201
202! **************************************************************************************************
203!> \brief ...
204! **************************************************************************************************
205 SUBROUTINE shuffle_check()
206 TYPE(rng_stream_type) :: rng_stream
207
208 INTEGER, PARAMETER :: sz = 20
209 INTEGER, DIMENSION(1:sz) :: arr, arr2, orig
210 LOGICAL, DIMENSION(1:sz) :: mask
211 INTEGER :: idx
212 REAL(KIND=dp), DIMENSION(3, 2), PARAMETER :: ig = 12.0_dp
213
214 WRITE (unit=default_output_unit, fmt="(/,/,T2,A)", advance="no") &
215 "Checking shuffle()"
216
217 rng_stream = rng_stream_type(name="shuffle() check", seed=ig)
218 orig = [(idx, idx=1, sz)]
219
220 arr = orig
221 CALL rng_stream%shuffle(arr)
222
223 IF (all(arr == orig)) then
224 error stop "shuffle failed: array was left untouched"
225 end if
226 WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
227
228 IF (any(arr /= orig(arr))) then
229 error stop "shuffle failed: the shuffled original is not the shuffled original"
230 end if
231 WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
232
233 ! sort and compare to orig
234 mask = .true.
235 DO idx = 1, size(orig)
236 IF (minval(arr, mask) /= orig(idx)) then
237 error stop "shuffle failed: there is at least one unknown index"
238 end if
239 mask(minloc(arr, mask)) = .false.
240 END DO
241 WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
242
243 arr2 = orig
244 CALL rng_stream%reset()
245 CALL rng_stream%shuffle(arr2)
246
247 IF (any(arr2 /= arr)) then
248 error stop "shuffle failed: array was shuffled differently with same rng state"
249 end if
250 WRITE (unit=default_output_unit, fmt="(A)", advance="no") "."
251
252 WRITE (unit=default_output_unit, fmt="(T4,A)") &
253 " successful"
254 END SUBROUTINE shuffle_check
255END PROGRAM parallel_rng_types_test
256! 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:46
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:141
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:25
message_passing::mp_world_init
subroutine, public mp_world_init(mp_comm)
initializes the system default communicator
Definition message_passing.F:1375
message_passing::mp_world_finalize
subroutine, public mp_world_finalize()
finalizes the system default communicator
Definition message_passing.F:1464
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:565
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:141
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:170
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157