gen_rand.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-10-18 15:09:55 pbrowne>
!!!
!!!    Collection of subroutines to make multidimensional random arrays
!!!    Copyright (C) 2014  Philip A. Browne
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program is distributed in the hope that it will be useful,
!!!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!!!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!!    GNU General Public License for more details.
!!!
!!!    You should have received a copy of the GNU General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: p.browne @ reading.ac.uk
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module random_number_controls
  character(10) :: normal_generator='random_d'
  contains
    subroutine set_random_number_controls
      use output_empire, only : emp_o,unit_nml,emp_e
      integer :: ios
      logical :: file_exists
      character(14) :: empire_namelist='empire.nml'
      
      namelist/random_number_controls/normal_generator
      
      inquire(file=empire_namelist,exist=file_exists)
      if(file_exists) then
         open(unit_nml,file=empire_namelist,iostat=ios,action='read'&
              &,status='old',form='formatted')
         if(ios .ne. 0) then
            write(emp_e,*) 'Cannot open ',empire_namelist
            write(emp_e,*) 'open_emp_o ERROR'
            stop
         end if
         read(unit_nml,nml=random_number_controls,iostat=ios)
         close(unit_nml)
         if(ios .ne. 0) normal_generator='random_d'
         write(emp_o,*) 'random_number_controls: normal_generator read as: &
              &',normal_generator
      else
         normal_generator = 'random_d'
      end if
    end subroutine set_random_number_controls
end module random_number_controls

Subroutine uniformrandomnumbers1d(minv, maxv, n,phi)
!use random
implicit none
integer, parameter :: rk = kind(1.0d0)
integer, intent(in) :: n
real(kind=rk), intent(in) :: minv
real(kind=rk), intent(in) :: maxv
real(kind=rk), dimension(n), intent(out) :: phi

call random_number(phi)

phi = minv + (maxv-minv)*phi
end Subroutine uniformrandomnumbers1d

Subroutine normalrandomnumbers1d(mean,stdev,n,phi)
  use output_empire, only : emp_e
use random_number_controls
use random
use ziggurat
IMPLICIT NONE
integer, parameter :: rk = kind(1.0d0)
integer, intent(in) :: n
real(kind=rk), INTENT(IN) :: mean
real(kind=rk), INTENT(IN) :: stdev
real(kind=rk), dimension(n), INTENT(OUT) :: phi
integer :: i

select case(normal_generator)
case('random_d')
   do i = 1,n
      phi(i) = mean+stdev*random_normal()
   end do
case('ziggurat')
   do i = 1,n
      phi(i) = mean+stdev*rnor()
   end do
case default
   write(emp_e,*) 'EMPIRE ERROR: wrong normal_generator selected in Nor&
        &malRandomNumbers1D'
   write(emp_e,*) 'EMPIRE ERROR: normal_generator = ',normal_generator&
        &,' is unknown'
   stop '-12'
end select



End Subroutine normalrandomnumbers1d

Subroutine normalrandomnumbers2d(mean,stdev,n,k,phi)
  use output_empire, only : emp_e
use random_number_controls
use ziggurat
use random
IMPLICIT NONE
integer, parameter :: rk = kind(1.0d0)
integer, intent(in) :: n
integer, intent(in) :: k
real(kind=rk), INTENT(IN) :: mean
real(kind=rk), INTENT(IN) :: stdev
real(kind=rk), dimension(n,k), INTENT(OUT) :: phi
integer :: i,j

select case(normal_generator)
case('random_d')
   do j = 1,k
      do i = 1,n
         phi(i,j) = mean+stdev*random_normal()
      end do
   end do
case('ziggurat')
   do j = 1,k
      do i = 1,n
         phi(i,j) = mean+stdev*rnor()
      end do
   end do
case default
   write(emp_e,*) 'EMPIRE ERROR: wrong normal_generator selected in Nor&
        &malRandomNumbers2D'
   write(emp_e,*) 'EMPIRE ERROR: normal_generator = ',normal_generator&
        &,' is unknown'
   stop '-12'
end select

End Subroutine normalrandomnumbers2d


subroutine mixturerandomnumbers1d(mean,stdev,ufac,epsi,n,phi,uniform)
use comms
use random
implicit none
include 'mpif.h'
real(kind=kind(1.0D0)), intent(in) :: mean,stdev,ufac,epsi
integer, intent(in) :: n
real(kind=kind(1.0D0)), dimension(n), intent(out) :: phi
logical, intent(out) :: uniform
real(kind=kind(1.0D0)) :: draw
integer :: mpi_err

if(comm_version .eq. 1 .or. comm_version .eq. 2) then
   call random_number(draw)
elseif(comm_version .eq. 3) then
   if(pf_member_rank .eq. 0) then
      call random_number(draw)
   end if
   call mpi_scatter(draw,1,mpi_double_precision,draw,1&
        &,mpi_double_precision,0,pf_member_rank,mpi_err)
else
    print*,'EMPIRE VERSION ',comm_version,' NOT SUPPORTED IN gen_rand'
    print*,'THIS IS AN ERROR. STOPPING'
    stop '-24'
end if

if(draw .gt. epsi) then
   call uniformrandomnumbers1d(mean-ufac,mean+ufac, n,phi)
   uniform = .true.
else
   call normalrandomnumbers1d(mean,stdev,n,phi)
   uniform = .false.
end if

end subroutine mixturerandomnumbers1d

subroutine mixturerandomnumbers2d(mean,stdev,ufac,epsi,n,k,phi,uniform)
use random
implicit none
real(kind=kind(1.0D0)), intent(in) :: mean,stdev,ufac,epsi
integer, intent(in) :: n,k
real(kind=kind(1.0D0)), dimension(n,k), intent(out) :: phi
logical, dimension(k), intent(out) :: uniform
real(kind=kind(1.0D0)) :: draw
integer :: i

do  i = 1,k
   call random_number(draw)

   if(draw .gt. epsi) then
      call uniformrandomnumbers1d(mean-ufac,mean+ufac,n,phi(:,i))
      uniform(i) = .true.
   else
      call normalrandomnumbers1d(mean,stdev,n,phi(:,i))
      uniform(i) = .false.
   end if
end do
end subroutine mixturerandomnumbers2d


subroutine random_seed_mpi(pfid)
  use random_number_controls
  use pf_control
  use ziggurat
  integer, intent(in) :: pfid

  integer :: n
  integer, allocatable, dimension(:) :: seed

  call set_random_number_controls

  call random_seed(size=n)
  allocate(seed(n))
  call random_seed(get=seed)
  !add the particle filter id to the seed to make it              
  !independent on each process              
  if(.not. pf%gen_data) then
     seed = seed+pfid
  else
     seed = seed-1
  end if
  call random_seed(put=seed)
  deallocate(seed)

  call zigset(pfid+1)

  

end subroutine random_seed_mpi