2 !-----------------------------------------------------------------------------
4 use prng ! prng.f90 or prng_32.f90
9 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
13 !-----------------------------------------------------------------------------
14 real(kind=8) function ran_number(x1,x2)
15 ! Calculate a random real number from the range (x1,x2).
17 use MPI_data ! include 'COMMON.SETUP'
21 ! implicit real*8 (a-h,o-z)
22 ! include 'DIMENSIONS'
23 real(kind=8) :: x1,x2,fctor
25 integer :: ix(1) !el ix ---> ix(1)
26 data fctor /2147483647.0D0/
28 ran_number=x1+(x2-x1)*prng_next(me)
31 ran_number=x1+(x2-x1)*ix(1)/fctor
34 end function ran_number
35 !-----------------------------------------------------------------------------
36 integer function iran_num(n1,n2)
37 ! Calculate a random integer number from the range (n1,n2).
39 use MPI_data ! include 'COMMON.SETUP'
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
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
53 ix(1)=n1+(n2-n1+1)*(ix(1)/fctor)
54 if (ix(1).gt.n2) ix(1)=n2
59 !-----------------------------------------------------------------------------
60 real(kind=8) function binorm(x1,x2,sigma1,sigma2,ak)
61 ! implicit real*8 (a-h,o-z)
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)
70 binorm=anorm_distr(x1,sigma1,alowb,aupb)
72 binorm=anorm_distr(x2,sigma2,alowb,aupb)
74 ! print '(a)','Exiting 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
83 real(kind=8) :: fac,gset,rsq,v1,v2,ran1
84 real(kind=8) :: x,sigma,alowb,aupb,gaussdev
87 !elwrite(iout,*) "anorm distr start",x,sigma,alowb,aupb
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
92 if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
93 fac=sqrt(-2.0d0*log(rsq)/rsq)
101 anorm_distr=x+gaussdev*sigma
102 !elwrite(iout,*) "anorm distr end",x,sigma,alowb,aupb,anorm_distr
104 end function anorm_distr
105 !-----------------------------------------------------------------------------
106 subroutine mult_norm(lda,n,a,x,fail)
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.
113 ! implicit double precision (a-h,o-z)
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
121 real(kind=8) :: sigma,alim,xi
123 ! print '(a)','Enter MULT_NORM.'
125 ! Find the smallest eigenvalue of the matrix A.
128 ! print '(8f10.5)',(a(i,j),j=1,n)
131 call f02faf('V','U',2,a,lda,eig,work,100,ifail)
133 call djacob(2,lda,10000,1.0d-10,a,vec,eig)
135 ! print '(8f10.5)',(eig(i),i=1,n)
138 ! print '(8f10.5)',(a(i,j),j=1,n)
140 if (eig(1).lt.eig_limit) then
141 print *,'From Mult_Norm: Eigenvalues of A are too small.'
146 ! Generate points following the normal distributions along the principal
147 ! axes of the moment matrix. Store in WORK.
150 sigma=1.0D0/dsqrt(eig(i))
152 work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
155 ! Transform the vector of normal variables back to the original basis.
165 end subroutine mult_norm
166 !-----------------------------------------------------------------------------
167 subroutine mult_norm1(lda,n,a,z,box,x,fail)
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.
175 ! implicit real*8 (a-h,o-z)
176 ! include 'DIMENSIONS'
178 use MPI_data !include 'COMMON.SETUP'
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'
191 ! Generate points following the normal distributions along the principal
192 ! axes of the moment matrix. Store in WORK.
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)
197 !d print *,i,box(1,i),box(2,i)
201 if (istep.gt.10000) then
202 ! write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
204 ! write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
206 ! write (iout,*) 'box',box
207 ! write (iout,*) 'a',a
208 ! write (iout,*) 'z',z
213 x(i)=ran_number(box(1,i),box(2,i))
218 ww=ww+0.5D0*a(i,i)*xi*xi
220 ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
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
231 if (dec.gt.etmp) goto 10
232 !d print *,'CG Processor',me,' exitting MultNorm1.'
234 end subroutine mult_norm1
235 !-----------------------------------------------------------------------------
237 !-----------------------------------------------------------------------------
240 integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
241 integer,parameter :: MASK=123459876
243 real(kind=4),parameter :: AM=1./IM
247 idum=IA*(idum-k*IQ)-IR*k
248 if (idum.lt.0) idum=idum+IM
253 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
254 !-----------------------------------------------------------------------------
257 integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
258 integer,parameter :: NTAB=32,NDIV=1+(IM-1)/NTAB
260 real(kind=4),parameter :: AM=1./IM,EPS=1.2e-7,RNMX=1.-EPS
262 integer,dimension(NTAB) :: iv
264 DATA iv /NTAB*0/, iy /0/
265 if (idum.le.0.or.iy.eq.0) then
269 idum=IA*(idum-k*IQ)-IR*k
270 if (idum.lt.0) idum=idum+IM
271 if (j.le.NTAB) iv(j)=idum
276 idum=IA*(idum-k*IQ)-IR*k
277 if (idum.lt.0) idum=idum+IM
284 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
285 !-----------------------------------------------------------------------------
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
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
296 DATA idum2/123456789/, iv/NTAB*0/, iy/0/
302 idum=IA1*(idum-k*IQ1)-k*IR1
303 if (idum.lt.0) idum=idum+IM1
304 if (j.le.NTAB) iv(j)=idum
309 idum=IA1*(idum-k*IQ1)-k*IR1
310 if (idum.lt.0) idum=idum+IM1
312 idum2=IA2*(idum2-k*IQ2)-k*IR2
313 if (idum2.lt.0) idum2=idum2+IM2
317 if(iy.lt.1)iy=iy+IMM1
321 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
322 !-----------------------------------------------------------------------------
325 integer,parameter :: MBIG=1000000000,MSEED=161803398,MZ=0
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
332 integer,dimension(55) :: ma
334 SAVE iff,inext,inextp,ma
336 if(idum.lt.0.or.iff.eq.0)then
346 if(mk.lt.MZ)mk=mk+MBIG
351 ma(i)=ma(i)-ma(1+mod(i+30,55))
352 if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
360 if(inext.eq.56)inext=1
362 if(inextp.eq.56)inextp=1
363 mj=ma(inext)-ma(inextp)
364 if(mj.lt.MZ)mj=mj+MBIG
369 ! (C) Copr. 1986-92 Numerical Recipes Software *11915
370 !-----------------------------------------------------------------------------
372 !-----------------------------------------------------------------------------
373 ! $Date: 1994/10/04 16:19:52 $
377 ! See help for RANDOMV on the PSFSHARE disk to understand these
378 ! subroutines. This is the VS Fortran version of this code.
381 subroutine VRND(VEC,N)
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
390 IF(.NOT.(I.GE.251))GOTO 23002
394 IF(.NOT.(I147.GE.251))GOTO 23004
397 A(I)=IEOR(A(I147),A(I))
402 !-----------------------------------------------------------------------------
403 real(kind=8) function RNDV(IDUM)
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
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
414 IX1 = MOD(IA1*IX1+IC1,M1)
416 IX1 = MOD(IA1*IX1+IC1,M1)
421 IX1 = MOD(IA1*IX1+IC1,M1)
422 IX2 = MOD(IA2*IX2+IC2,M2)
423 R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
426 IX1 = MOD(IA1*IX1+IC1,M1)
427 IX2 = MOD(IA2*IX2+IC2,M2)
428 IX3 = MOD(IA3*IX3+IC3,M3)
431 R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
435 !-----------------------------------------------------------------------------
436 subroutine VRNDST(SEED)
439 !el INTEGER :: A(250)
440 INTEGER :: LOOP,IDUM !el,I,I147
442 !el real(kind=8) :: RNDV
443 !el COMMON /VRANDD/ A, I, I147
448 A(LOOP)=INT(RNDV(IDUM)*2147483647)
451 end subroutine VRNDST
452 !-----------------------------------------------------------------------------
453 subroutine VRNDIN(IODEV)
455 INTEGER :: IODEV !el, A(250),I, I147
456 !el COMMON/VRANDD/ A, I, I147
457 READ(IODEV) A, I, I147
459 end subroutine VRNDIN
460 !-----------------------------------------------------------------------------
461 subroutine VRNDOU(IODEV)
462 ! This corresponds to VRNDOUT in the APFTN64 version
464 INTEGER :: IODEV !el, A(250),I, I147
465 !el COMMON/VRANDD/ A, I, I147
466 WRITE(IODEV) A, I, I147
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
476 RNUNF= DBLE( IRAN1(1) ) / FCTOR
477 !******************************
478 ! write(6,*) 'rnunf in rnunf = ',rnunf
481 !-----------------------------------------------------------------------------
483 !-----------------------------------------------------------------------------
484 subroutine random_init(seed)
486 ! Initialize random number generator
488 ! implicit real*8 (a-h,o-z)
489 ! include 'DIMENSIONS'
493 logical :: OKRandom !, prng_restart
495 integer :: iseed_array(4)
496 integer(kind=8) :: iseed
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))
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.'
523 call mpi_finalize(mpi_comm_world,ierr)
525 stop 'Bad random seed.'
528 if (fg_rank.eq.0) then
532 write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
533 OKRandom = prng_restart(me,iseed)
537 iseed_array(i) = dint(seed/tmp)
538 seed=seed-iseed_array(i)*tmp
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)
549 r1=ran_number(0.0D0,1.0D0)
551 write (iout,*) 'ran_num',r1
552 if (r1.lt.0.0d0) OKRandom=.false.
554 if (.not.OKRandom) then
555 write (iout,*) 'PRNG IS NOT WORKING!!!'
556 print *,'PRNG IS NOT WORKING!!!'
559 call mpi_abort(mpi_comm_world,error_msg,ierr)
562 write (iout,*) 'too many processors for parallel prng'
563 write (*,*) 'too many processors for parallel prng'
571 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
574 end subroutine random_init
575 !-----------------------------------------------------------------------------
576 !-----------------------------------------------------------------------------