./28_numeric_noabirule/uniformrandom.F90

ABINIT/uniformrandom

ABINIT/uniformrandom [ Functions ]

NAME

  uniformrandom

FUNCTION

  Copyright (C) 2014-2017 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

 28 #if defined HAVE_CONFIG_H
 29 #include "config.h"
 30 #endif
 31 
 32 #include "abi_common.h"
 33 
 34 function uniformrandom(seed) 
 35 ! Returns a uniform random deviate between 0.0 and 1.0.
 36 ! Set seed to any value < 0 to initialize or reinitialize sequence.
 37 ! Parameters are chosen from integer overflow=2**23 (conservative).
 38 ! For some documentation, see Numerical Recipes, 1986, p196.
 39  use m_errors
 40 
 41 !This section has been created automatically by the script Abilint (TD).
 42 !Do not modify the following lines by hand.
 43 #undef ABI_FUNC
 44 #define ABI_FUNC 'uniformrandom'
 45 !End of the abilint section
 46 
 47  implicit none
 48 !Input/Output
 49  double precision :: uniformrandom
 50  integer :: seed
 51 !Local variables
 52  integer, parameter :: im1=11979,ia1= 430,ic1=2531
 53  integer, parameter :: im2= 6655,ia2= 936,ic2=1399
 54  integer, parameter :: im3= 6075,ia3=1366,ic3=1283
 55  integer, save :: init=0
 56  integer, save :: ii1,ii2,ii3
 57  integer :: kk 
 58  double precision :: im1inv,im2inv
 59  double precision, save :: table(97)
 60  character(len=500) :: message
 61 
 62  im1inv=1.0d0/im1 ; im2inv=1.0d0/im2
 63 
 64 !Initialize on first call or when seed<0:
 65  if (seed<0.or.init==0) then
 66    seed=-abs(seed) 
 67 
 68 !  First generator
 69    ii1=mod(ic1-seed,im1)
 70    ii1=mod(ia1*ii1+ic1,im1) 
 71 !  Second generator
 72    ii2=mod(ii1,im2)
 73    ii1=mod(ia1*ii1+ic1,im1) 
 74 !  Third generator
 75    ii3=mod(ii1,im3) 
 76 
 77 !  Fill table
 78    do kk=1,97
 79      ii1=mod(ia1*ii1+ic1,im1)
 80      ii2=mod(ia2*ii2+ic2,im2) 
 81      table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
 82    enddo
 83 
 84    init=1 ; seed=1
 85  end if 
 86 
 87 !Third generator gives index
 88  ii3=mod(ia3*ii3+ic3,im3) 
 89  kk=1+(97*ii3)/im3
 90  if (kk<1.or.kk>97) then
 91    write(message,'(a,2i0,a)' ) ' trouble in uniformrandom; ii3,kk=',ii3,kk,' =>stop'
 92    MSG_ERROR(message)
 93  end if 
 94  uniformrandom=table(kk) 
 95 
 96 !Replace old value, based on generators 1 and 2
 97  ii1=mod(ia1*ii1+ic1,im1)
 98  ii2=mod(ia2*ii2+ic2,im2)
 99  table(kk)=(dble(ii1)+dble(ii2)*im2inv)*im1inv
100 
101 end function uniformrandom

TABLE OF CONTENTS

ABINIT/uniformrandom [ Functions ]