source/unres/random.f90

   1       module random
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use prng  ! prng.f90 or prng_32.f90
   5       use math
   6       implicit none
   7 !
   8 !-----------------------------------------------------------------------------
   9       contains
  10 !-----------------------------------------------------------------------------
  11 ! gen_rand_conf.F
  12 !-----------------------------------------------------------------------------
  13       real(kind=8) function ran_number(x1,x2)
  14 ! Calculate a random real number from the range (x1,x2).
  15 #ifdef MPI
  16       use MPI_data      ! include 'COMMON.SETUP'
  17       include "mpif.h"
  18 #endif
  19
  20 !      implicit real*8 (a-h,o-z)
  21 !      include 'DIMENSIONS'
  22       real(kind=8) :: x1,x2,fctor
  23 !el local variables
  24       integer :: ix(1)  !el ix --->  ix(1)
  25       data fctor /2147483647.0D0/
  26 #ifdef MPI
  27       ran_number=x1+(x2-x1)*prng_next(me)
  28 #else
  29       call vrnd(ix(1),1)
  30       ran_number=x1+(x2-x1)*ix(1)/fctor
  31 #endif
  32       return
  33       end function ran_number
  34 !-----------------------------------------------------------------------------
  35       integer function iran_num(n1,n2)
  36 ! Calculate a random integer number from the range (n1,n2).
  37 #ifdef MPI
  38       use MPI_data      ! include 'COMMON.SETUP'
  39       include "mpif.h"
  40 #endif
  41 !      implicit real*8 (a-h,o-z)
  42 !      include 'DIMENSIONS'
  43       integer :: n1,n2,ix(1)    !el ix --->  ix(1)
  44       real(kind=4) :: fctor=2147483647.0
  45 #ifdef MPI
  46       ix(1)=n1+(n2-n1+1)*prng_next(me)
  47       if (ix(1).lt.n1) ix(1)=n1
  48       if (ix(1).gt.n2) ix(1)=n2
  49       iran_num=ix(1)
  50 #else
  51       call vrnd(ix(1),1)
  52       ix(1)=n1+(n2-n1+1)*(ix(1)/fctor)
  53       if (ix(1).gt.n2) ix(1)=n2
  54       iran_num=ix(1)
  55 #endif
  56       return
  57       end function iran_num
  58 !-----------------------------------------------------------------------------
  59       real(kind=8) function binorm(x1,x2,sigma1,sigma2,ak)
  60 !      implicit real*8 (a-h,o-z)
  61 !el local variables
  62       real(kind=8) :: x1,x2,sigma1,sigma2,ak,alowb,aupb,seg,alen
  63 !     print '(a)','Enter BINORM.'
  64       alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
  65       aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
  66       seg=sigma1/(sigma1+ak*sigma2)
  67       alen=ran_number(0.0D0,1.0D0)
  68       if (alen.lt.seg) then
  69         binorm=anorm_distr(x1,sigma1,alowb,aupb)
  70       else
  71         binorm=anorm_distr(x2,sigma2,alowb,aupb)
  72       endif
  73 !     print '(a)','Exiting BINORM.'
  74       return
  75       end function binorm
  76 !-----------------------------------------------------------------------------
  77       real(kind=8) function anorm_distr(x,sigma,alowb,aupb)
  78 !        implicit real*8 (a-h,o-z)
  79 !  to make a normally distributed deviate with zero mean and unit variance
  80 !
  81         integer :: iset
  82         real(kind=8) :: fac,gset,rsq,v1,v2,ran1
  83         real(kind=8) :: x,sigma,alowb,aupb,gaussdev
  84         save iset,gset
  85         data iset/0/
  86 !elwrite(iout,*) "anorm distr start",x,sigma,alowb,aupb
  87         if(iset.eq.0) then
  88 1               v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
  89                 v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
  90                 rsq=v1**2+v2**2
  91                 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
  92                 fac=sqrt(-2.0d0*log(rsq)/rsq)
  93                 gset=v1*fac
  94                 gaussdev=v2*fac
  95                 iset=1
  96         else
  97                 gaussdev=gset
  98                 iset=0
  99         endif
 100         anorm_distr=x+gaussdev*sigma
 101 !elwrite(iout,*) "anorm distr end",x,sigma,alowb,aupb,anorm_distr
 102       return
 103       end function anorm_distr
 104 !-----------------------------------------------------------------------------
 105       subroutine mult_norm(lda,n,a,x,fail)
 106 !
 107 ! Generate the vector X whose elements obey the multiple-normal distribution
 108 ! from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
 109 ! n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
 110 ! eigenvalue of the matrix A is close to 0.
 111 !
 112 !      implicit double precision (a-h,o-z)
 113       integer :: lda,n,i,j
 114       real(kind=8),dimension(lda,n) :: a
 115       real(kind=8),dimension(n) :: x
 116       real(kind=8),dimension(100) :: eig,work
 117       real(kind=8),dimension(3,3) :: vec
 118       real(kind=8) :: eig_limit=1.0D-8
 119       logical :: fail
 120       real(kind=8) :: sigma,alim,xi
 121       fail=.false.
 122 !     print '(a)','Enter MULT_NORM.'
 123 !
 124 ! Find the smallest eigenvalue of the matrix A.
 125 !
 126 !     do i=1,n
 127 !       print '(8f10.5)',(a(i,j),j=1,n)
 128 !     enddo
 129 #ifdef NAG
 130       call f02faf('V','U',2,a,lda,eig,work,100,ifail)
 131 #else
 132       call djacob(2,lda,10000,1.0d-10,a,vec,eig)
 133 #endif
 134 !     print '(8f10.5)',(eig(i),i=1,n)
 135 !     print '(a)'
 136 !     do i=1,n
 137 !       print '(8f10.5)',(a(i,j),j=1,n)
 138 !     enddo
 139       if (eig(1).lt.eig_limit) then
 140         print *,'From Mult_Norm: Eigenvalues of A are too small.'
 141         fail=.true.
 142         return
 143       endif
 144 !
 145 ! Generate points following the normal distributions along the principal
 146 ! axes of the moment matrix. Store in WORK.
 147 !
 148       do i=1,n
 149         sigma=1.0D0/dsqrt(eig(i))
 150         alim=-3.0D0*sigma
 151         work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
 152       enddo
 153 !
 154 ! Transform the vector of normal variables back to the original basis.
 155 !
 156       do i=1,n
 157         xi=0.0D0
 158         do j=1,n
 159           xi=xi+a(i,j)*work(j)
 160         enddo
 161         x(i)=xi
 162       enddo
 163       return
 164       end subroutine mult_norm
 165 !-----------------------------------------------------------------------------
 166       subroutine mult_norm1(lda,n,a,z,box,x,fail)
 167 !
 168 ! Generate the vector X whose elements obey the multi-gaussian multi-dimensional
 169 ! distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
 170 ! leading dimension of the moment matrix A, n is the dimension of the
 171 ! distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
 172 ! smallest eigenvalue of the matrix A is close to 0.
 173 !
 174 !      implicit real*8 (a-h,o-z)
 175 !      include 'DIMENSIONS'
 176 #ifdef MP
 177       use MPI_data      !include 'COMMON.SETUP'
 178 #endif
 179 #ifdef MPI
 180       include 'mpif.h'
 181 #endif
 182       integer :: lda,n,i,j,istep
 183       real(kind=8),dimension(lda,n) :: a
 184       real(kind=8),dimension(n) :: z,x
 185       real(kind=8),dimension(n,n) :: box
 186       real(kind=8) :: etmp,ww,xi,dec
 187 !      include 'COMMON.IOUNITS'
 188       logical :: fail
 189 !
 190 ! Generate points following the normal distributions along the principal
 191 ! axes of the moment matrix. Store in WORK.
 192 !
 193 !d    print *,'CG Processor',me,' entered MultNorm1.'
 194 !d    print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
 195 !d    do i=1,n
 196 !d      print *,i,box(1,i),box(2,i)
 197 !d    enddo
 198       istep = 0
 199    10 istep = istep + 1
 200       if (istep.gt.10000) then
 201 !        write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 202 !     & ' in MultNorm1.'
 203 !        write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
 204 !     & ' in MultNorm1.'
 205 !        write (iout,*) 'box',box
 206 !        write (iout,*) 'a',a
 207 !        write (iout,*) 'z',z
 208         fail=.true.
 209         return
 210       endif
 211       do i=1,n
 212        x(i)=ran_number(box(1,i),box(2,i))
 213       enddo
 214       ww=0.0D0
 215       do i=1,n
 216         xi=pinorm(x(i)-z(i))
 217         ww=ww+0.5D0*a(i,i)*xi*xi
 218         do j=i+1,n
 219           ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
 220         enddo
 221       enddo
 222       dec=ran_number(0.0D0,1.0D0)
 223 !      print *,(x(i),i=1,n),ww,dexp(-ww),dec
 224 !rc   if (dec.gt.dexp(-ww)) goto 10
 225       if(-ww.lt.100) then
 226        etmp=dexp(-ww)
 227       else
 228        return
 229       endif
 230       if (dec.gt.etmp) goto 10
 231 !d    print *,'CG Processor',me,' exitting MultNorm1.'
 232       return
 233       end subroutine mult_norm1
 234 !-----------------------------------------------------------------------------
 235 ! ran.f
 236 !-----------------------------------------------------------------------------
 237       function ran0(idum)
 238       integer :: idum
 239       integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
 240       integer,parameter :: MASK=123459876
 241       real(kind=4) :: ran0
 242       real(kind=4),parameter :: AM=1./IM
 243       integer :: k
 244       idum=ieor(idum,MASK)
 245       k=idum/IQ
 246       idum=IA*(idum-k*IQ)-IR*k
 247       if (idum.lt.0) idum=idum+IM
 248       ran0=AM*idum
 249       idum=ieor(idum,MASK)
 250       return
 251       end function ran0
 252 !  (C) Copr. 1986-92 Numerical Recipes Software *11915
 253 !-----------------------------------------------------------------------------
 254       function ran1(idum)
 255       integer :: idum
 256       integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
 257       integer,parameter :: NTAB=32,NDIV=1+(IM-1)/NTAB
 258       real(kind=4) :: ran1
 259       real(kind=4),parameter :: AM=1./IM,EPS=1.2e-7,RNMX=1.-EPS
 260       integer :: j,k,iy
 261       integer,dimension(NTAB) :: iv
 262       SAVE iv,iy
 263       DATA iv /NTAB*0/, iy /0/
 264       if (idum.le.0.or.iy.eq.0) then
 265         idum=max(-idum,1)
 266         do 11 j=NTAB+8,1,-1
 267           k=idum/IQ
 268           idum=IA*(idum-k*IQ)-IR*k
 269           if (idum.lt.0) idum=idum+IM
 270           if (j.le.NTAB) iv(j)=idum
 271 11      continue
 272         iy=iv(1)
 273       endif
 274       k=idum/IQ
 275       idum=IA*(idum-k*IQ)-IR*k
 276       if (idum.lt.0) idum=idum+IM
 277       j=1+iy/NDIV
 278       iy=iv(j)
 279       iv(j)=idum
 280       ran1=min(AM*iy,RNMX)
 281       return
 282       end function ran1
 283 !  (C) Copr. 1986-92 Numerical Recipes Software *11915
 284 !-----------------------------------------------------------------------------
 285       function ran2(idum)
 286       integer :: idum
 287       integer,parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1,&
 288               IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
 289               NTAB=32,NDIV=1+IMM1/NTAB
 290       real(kind=4) :: ran2
 291       real(kind=4),parameter :: AM=1./IM1,EPS=1.2e-7,RNMX=1.-EPS
 292       integer :: idum2,j,k,iy
 293       integer,dimension(NTAB) :: iv
 294       SAVE iv,iy,idum2
 295       DATA idum2/123456789/, iv/NTAB*0/, iy/0/
 296       if (idum.le.0) then
 297         idum=max(-idum,1)
 298         idum2=idum
 299         do 11 j=NTAB+8,1,-1
 300           k=idum/IQ1
 301           idum=IA1*(idum-k*IQ1)-k*IR1
 302           if (idum.lt.0) idum=idum+IM1
 303           if (j.le.NTAB) iv(j)=idum
 304 11      continue
 305         iy=iv(1)
 306       endif
 307       k=idum/IQ1
 308       idum=IA1*(idum-k*IQ1)-k*IR1
 309       if (idum.lt.0) idum=idum+IM1
 310       k=idum2/IQ2
 311       idum2=IA2*(idum2-k*IQ2)-k*IR2
 312       if (idum2.lt.0) idum2=idum2+IM2
 313       j=1+iy/NDIV
 314       iy=iv(j)-idum2
 315       iv(j)=idum
 316       if(iy.lt.1)iy=iy+IMM1
 317       ran2=min(AM*iy,RNMX)
 318       return
 319       end function ran2
 320 !  (C) Copr. 1986-92 Numerical Recipes Software *11915
 321 !-----------------------------------------------------------------------------
 322       function ran3(idum)
 323       integer :: idum
 324       integer,parameter :: MBIG=1000000000,MSEED=161803398,MZ=0
 325 !     REAL MBIG,MSEED,MZ
 326       real(kind=4) :: ran3
 327       real(kind=4),parameter :: FAC=1./MBIG
 328 !     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
 329       integer :: i,iff,ii,inext,inextp,k
 330       integer :: mj,mk
 331       integer,dimension(55) :: ma
 332 !     REAL mj,mk,ma(55)
 333       SAVE iff,inext,inextp,ma
 334       DATA iff /0/
 335       if(idum.lt.0.or.iff.eq.0)then
 336         iff=1
 337         mj=MSEED-iabs(idum)
 338         mj=mod(mj,MBIG)
 339         ma(55)=mj
 340         mk=1
 341         do 11 i=1,54
 342           ii=mod(21*i,55)
 343           ma(ii)=mk
 344           mk=mj-mk
 345           if(mk.lt.MZ)mk=mk+MBIG
 346           mj=ma(ii)
 347 11      continue
 348         do 13 k=1,4
 349           do 12 i=1,55
 350             ma(i)=ma(i)-ma(1+mod(i+30,55))
 351             if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
 352 12        continue
 353 13      continue
 354         inext=0
 355         inextp=31
 356         idum=1
 357       endif
 358       inext=inext+1
 359       if(inext.eq.56)inext=1
 360       inextp=inextp+1
 361       if(inextp.eq.56)inextp=1
 362       mj=ma(inext)-ma(inextp)
 363       if(mj.lt.MZ)mj=mj+MBIG
 364       ma(inext)=mj
 365       ran3=mj*FAC
 366       return
 367       end function ran3
 368 !  (C) Copr. 1986-92 Numerical Recipes Software *11915
 369 !-----------------------------------------------------------------------------
 370 ! randgens.f
 371 !-----------------------------------------------------------------------------
 372 ! $Date: 1994/10/04 16:19:52 $
 373 ! $Revision: 2.1 $
 374 !
 375 !
 376 !  See help for RANDOMV on the PSFSHARE disk to understand these
 377 !  subroutines.  This is the VS Fortran version of this code.
 378 !
 379 !
 380       subroutine VRND(VEC,N)
 381
 382       use comm_vrandd
 383 !el      integer,dimension(250) :: A
 384       integer :: LOOP,N !el,I,I147
 385       integer,dimension(N) :: VEC
 386 !el      COMMON /VRANDD/ A, I, I147
 387       DO 23000 LOOP=1,N
 388       I=I+1
 389       IF(.NOT.(I.GE.251))GOTO 23002
 390       I=1
 391 23002 CONTINUE
 392       I147=I147+1
 393       IF(.NOT.(I147.GE.251))GOTO 23004
 394       I147=1
 395 23004 CONTINUE
 396       A(I)=IEOR(A(I147),A(I))
 397       VEC(LOOP)=A(I)
 398 23000 CONTINUE
 399       return
 400       end subroutine VRND
 401 !-----------------------------------------------------------------------------
 402       real(kind=8) function RNDV(IDUM)
 403
 404       real(kind=8) :: RM1,RM2,R(99)
 405       INTEGER :: IA1,IC1,M1, IA2,IC2,M2, IA3,IC3,M3, IDUM
 406       INTEGER :: IX1,IX2,IX3,J
 407       SAVE
 408       DATA IA1,IC1,M1/1279,351762,1664557/
 409       DATA IA2,IC2,M2/2011,221592,1048583/
 410       DATA IA3,IC3,M3/15551,6150,29101/
 411       IF(.NOT.(IDUM.LT.0))GOTO 23006
 412       IX1 = MOD(-IDUM,M1)
 413       IX1 = MOD(IA1*IX1+IC1,M1)
 414       IX2 = MOD(IX1,M2)
 415       IX1 = MOD(IA1*IX1+IC1,M1)
 416       IX3 = MOD(IX1,M3)
 417       RM1 = 1./DBLE(M1)
 418       RM2 = 1./DBLE(M2)
 419       DO 23008 J = 1,99
 420       IX1 = MOD(IA1*IX1+IC1,M1)
 421       IX2 = MOD(IA2*IX2+IC2,M2)
 422       R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
 423 23008 CONTINUE
 424 23006 CONTINUE
 425       IX1 = MOD(IA1*IX1+IC1,M1)
 426       IX2 = MOD(IA2*IX2+IC2,M2)
 427       IX3 = MOD(IA3*IX3+IC3,M3)
 428       J = 1+(99*IX3)/M3
 429       RNDV = R(J)
 430       R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
 431       IDUM = IX1
 432       return
 433       end function RNDV
 434 !-----------------------------------------------------------------------------
 435       subroutine VRNDST(SEED)
 436
 437       use comm_vrandd
 438 !el      INTEGER :: A(250)
 439       INTEGER :: LOOP,IDUM !el,I,I147
 440       INTEGER :: SEED
 441 !el      real(kind=8) :: RNDV
 442 !el      COMMON /VRANDD/ A, I, I147
 443       I=0
 444       I147=103
 445       IDUM=SEED
 446       DO 23010 LOOP=1,250
 447       A(LOOP)=INT(RNDV(IDUM)*2147483647)
 448 23010 CONTINUE
 449       return
 450       end subroutine VRNDST
 451 !-----------------------------------------------------------------------------
 452       subroutine VRNDIN(IODEV)
 453       use comm_vrandd
 454       INTEGER :: IODEV !el, A(250),I, I147
 455 !el      COMMON/VRANDD/ A, I, I147
 456       READ(IODEV) A, I, I147
 457       return
 458       end subroutine VRNDIN
 459 !-----------------------------------------------------------------------------
 460       subroutine VRNDOU(IODEV)
 461 !       This corresponds to VRNDOUT in the APFTN64 version
 462       use comm_vrandd
 463       INTEGER :: IODEV !el, A(250),I, I147
 464 !el      COMMON/VRANDD/ A, I, I147
 465       WRITE(IODEV) A, I, I147
 466       return
 467       end subroutine VRNDOU
 468 !-----------------------------------------------------------------------------
 469       real(kind=8) function RNUNF(N)
 470       INTEGER :: IRAN1(2000),N
 471       real(kind=8) :: FCTOR
 472       DATA FCTOR /2147483647.0D0/
 473 !     We get only one random number, here!    DR  9/1/92
 474       CALL VRND(IRAN1,1)
 475       RNUNF= DBLE( IRAN1(1) ) / FCTOR
 476 !******************************
 477 !     write(6,*) 'rnunf  in rnunf = ',rnunf
 478       return
 479       end function RNUNF
 480 !-----------------------------------------------------------------------------
 481 ! readrtns_CSA.F
 482 !-----------------------------------------------------------------------------
 483       subroutine random_init(seed)
 484 !
 485 ! Initialize random number generator
 486 !
 487 !      implicit real*8 (a-h,o-z)
 488 !      include 'DIMENSIONS'
 489 #ifdef MPI
 490       use MPI_data
 491       include 'mpif.h'
 492       logical :: OKRandom !, prng_restart
 493       real(kind=8) :: r1
 494       integer :: iseed_array(4)
 495       integer(kind=8) :: iseed
 496 #else
 497       integer :: iseed
 498 #endif
 499       real(kind=8) :: seed
 500       integer :: ierr,error_msg
 501 !      include 'COMMON.IOUNITS'
 502 !      include 'COMMON.TIME1'
 503 !      include 'COMMON.THREAD'
 504 !      include 'COMMON.SBRIDGE'
 505 !      include 'COMMON.CONTROL'
 506 !      include 'COMMON.MCM'
 507 !      include 'COMMON.MAP'
 508 !      include 'COMMON.HEADER'
 509 !      include 'COMMON.CSA'
 510 !      include 'COMMON.CHAIN'
 511 !      include 'COMMON.MUCA'
 512 !      include 'COMMON.MD'
 513 !      include 'COMMON.FFIELD'
 514 !      include 'COMMON.SETUP'
 515       iseed=-dint(dabs(seed))
 516       if (iseed.eq.0) then
 517         write (iout,'(/80(1h*)/20x,a/80(1h*))') &
 518           'Random seed undefined. The program will stop.'
 519         write (*,'(/80(1h*)/20x,a/80(1h*))') &
 520           'Random seed undefined. The program will stop.'
 521 #ifdef MPI
 522         call mpi_finalize(mpi_comm_world,ierr)
 523 #endif
 524         stop 'Bad random seed.'
 525       endif
 526 #ifdef MPI
 527       if (fg_rank.eq.0) then
 528       seed=seed*(me+1)+1
 529 #ifdef AMD64
 530       if(me.eq.king) &
 531         write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
 532       OKRandom = prng_restart(me,iseed)
 533 #else
 534       do i=1,4
 535        tmp=65536.0d0**(4-i)
 536        iseed_array(i) = dint(seed/tmp)
 537        seed=seed-iseed_array(i)*tmp
 538       enddo
 539       if(me.eq.king) &
 540        write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',&
 541                        (iseed_array(i),i=1,4)
 542       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',&
 543                        (iseed_array(i),i=1,4)
 544       OKRandom = prng_restart(me,iseed_array)
 545 #endif
 546       if (OKRandom) then
 547 !        r1 = prng_next(me)
 548          r1=ran_number(0.0D0,1.0D0)
 549         if(me.eq.king) &
 550          write (iout,*) 'ran_num',r1
 551         if (r1.lt.0.0d0) OKRandom=.false.
 552       endif
 553       if (.not.OKRandom) then
 554         write (iout,*) 'PRNG IS NOT WORKING!!!'
 555         print *,'PRNG IS NOT WORKING!!!'
 556         if (me.eq.0) then
 557          call flush(iout)
 558          call mpi_abort(mpi_comm_world,error_msg,ierr)
 559          stop
 560         else
 561          write (iout,*) 'too many processors for parallel prng'
 562          write (*,*) 'too many processors for parallel prng'
 563          call flush(iout)
 564          stop
 565         endif
 566       endif
 567       endif
 568 #else
 569       call vrndst(iseed)
 570       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
 571 #endif
 572       return
 573       end subroutine random_init
 574 !-----------------------------------------------------------------------------
 575 !-----------------------------------------------------------------------------
 576       end module random