DISPATCH
random_mod.f90
1 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2 !> $Id: 2219c9172a3fc9232139b1fbb4bed9faac2e5d13 $
3 !> Simple random number module
4 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
5 MODULE random_mod
6  USE io_mod
7  USE trace_mod
8  implicit none
9  private
10  interface uniform
11  module procedure uniform_1, uniform_n
12  end interface
13  type, public:: random_t
14  integer:: seed=11
15  integer:: id
16  contains
17  procedure:: init
18  procedure:: ran1
19  procedure:: ran3
20  procedure:: uniform_1
21  procedure:: uniform_n
22  end type
23  integer, save:: id=0
24 CONTAINS
25 
26 !===============================================================================
27 !> Initialize both random number generators
28 !>
29 !> Note that, if the random number is use to create a coherent pattern as a
30 !> function of code time, then -- since the differen tasks are in general not at
31 !> the same code time, each task needs to use the task%random instance, with
32 !> the same seed, but updated independently, as a function of their local code
33 !> time. For this reason, there is no generic, shared instance of the random_t
34 !> data type.
35 !===============================================================================
36 SUBROUTINE init (self, seed)
37  class(random_t):: self
38  integer, optional:: seed
39  real:: a, b(3)
40  !.............................................................................
41  call trace_begin ('random_t%init')
42  !$omp critical (random_cr)
43  id = id+1
44  self%id = id
45  !$omp end critical (random_cr)
46  if (present(seed)) then
47  self%seed = abs(seed)
48  else
49  self%seed = 11
50  end if
51  call random_seed
52  !call random_number (a)
53  !call random_number (b)
54  !if (io%master) then
55  ! print*,'random_number:',a
56  ! print*,'random_number:',b
57  !end if
58  call trace_end
59 END SUBROUTINE init
60 
61 !===============================================================================
62 !> Ancient Numerical Recipes code. Will do for now.
63 !===============================================================================
64 FUNCTION ran1(self)
65  class(random_t):: self
66  integer:: idum,ia,im,iq,ir,ntab,ndiv
67  real:: ran1,am,eps,rnmx
68  parameter(ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836, &
69  ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
70  integer:: k,iy
71  !.............................................................................
72  if (self%id==0) &
73  print *,'WARNING: ran1 has not been initialized'
74  idum = self%seed
75  k =idum/iq
76  idum = ia*(idum-k*iq)-ir*k
77  if (idum<0) idum = idum+im
78  iy = idum
79  ran1 = min(am*iy,rnmx)
80  self%seed = idum
81 END FUNCTION
82 
83 !===============================================================================
84 !> 3-component random vector, components uniform in a spherical distribution,
85 !> with max radius = 1.0
86 !===============================================================================
87 FUNCTION ran3(self) result (u)
88  class(random_t):: self
89  real(kind=4):: u(3), u2
90  integer:: i
91  !-----------------------------------------------------------------------------
92  do i=1,3
93  u(i) = 2.*self%ran1()-1.
94  end do
95  u2 = sum(u**2)
96  do while (u2 > 1.0 .or. u2 < 1e-8)
97  do i=1,3
98  u(i) = 2.*self%ran1()-1.
99  end do
100  u2 = sum(u**2)
101  end do
102 END FUNCTION ran3
103 
104 !===============================================================================
105 !===============================================================================
106 FUNCTION uniform_1 (self, a, b) RESULT (out)
107  class(random_t):: self
108  real, intent(in), optional:: a, b
109  real:: out
110  !-----------------------------------------------------------------------------
111  if (present(a)) then
112  if (present(b)) then
113  out = a + (b-a)*self%ran1()
114  else
115  out = a*self%ran1()
116  end if
117  else
118  out = self%ran1()
119  end if
120 END FUNCTION
121 
122 !===============================================================================
123 !===============================================================================
124 FUNCTION uniform_n (self, n, a, b) RESULT (out)
125  class(random_t):: self
126  integer, intent(in):: n
127  real, intent(in), optional:: a, b
128  real, dimension(n):: out
129  integer:: i
130  !-----------------------------------------------------------------------------
131  if (present(a)) then
132  if (present(b)) then
133  do i=1,n
134  out(i) = a*self%ran1()
135  end do
136  else
137  do i=1,n
138  out(i) = a + (b-a)*self%ran1()
139  end do
140  end if
141  else
142  do i=1,n
143  out(i) = self%ran1()
144  end do
145  end if
146 END FUNCTION
147 
148 END MODULE random_mod
Definition: io_mod.f90:4