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