+++ /dev/null
- module random
-!-----------------------------------------------------------------------------
- use io_units
- use prng ! prng.f90 or prng_32.f90
- use math
- implicit none
-! public :: rndv
-!
-!-----------------------------------------------------------------------------
- contains
-!-----------------------------------------------------------------------------
-! gen_rand_conf.F
-!-----------------------------------------------------------------------------
- real(kind=8) function ran_number(x1,x2)
-! Calculate a random real number from the range (x1,x2).
-#ifdef MPI
- use MPI_data ! include 'COMMON.SETUP'
- include "mpif.h"
-#endif
-
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- real(kind=8) :: x1,x2,fctor
-!el local variables
- integer :: ix(1) !el ix ---> ix(1)
- data fctor /2147483647.0D0/
-#ifdef MPI
- ran_number=x1+(x2-x1)*prng_next(me)
-#else
- call vrnd(ix(1),1)
- ran_number=x1+(x2-x1)*ix(1)/fctor
-#endif
- return
- end function ran_number
-!-----------------------------------------------------------------------------
- integer function iran_num(n1,n2)
-! Calculate a random integer number from the range (n1,n2).
-#ifdef MPI
- use MPI_data ! include 'COMMON.SETUP'
- include "mpif.h"
-#endif
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- integer :: n1,n2,ix(1) !el ix ---> ix(1)
- real(kind=4) :: fctor=2147483647.0
-#ifdef MPI
- ix(1)=n1+(n2-n1+1)*prng_next(me)
- if (ix(1).lt.n1) ix(1)=n1
- if (ix(1).gt.n2) ix(1)=n2
- iran_num=ix(1)
-#else
- call vrnd(ix(1),1)
- ix(1)=n1+(n2-n1+1)*(ix(1)/fctor)
- if (ix(1).gt.n2) ix(1)=n2
- iran_num=ix(1)
-#endif
- return
- end function iran_num
-!-----------------------------------------------------------------------------
- real(kind=8) function binorm(x1,x2,sigma1,sigma2,ak)
-! implicit real*8 (a-h,o-z)
-!el local variables
- real(kind=8) :: x1,x2,sigma1,sigma2,ak,alowb,aupb,seg,alen
-! print '(a)','Enter BINORM.'
- alowb=dmin1(x1-3.0D0*sigma1,x2-3.0D0*sigma2)
- aupb=dmax1(x1+3.0D0*sigma1,x2+3.0D0*sigma2)
- seg=sigma1/(sigma1+ak*sigma2)
- alen=ran_number(0.0D0,1.0D0)
- if (alen.lt.seg) then
- binorm=anorm_distr(x1,sigma1,alowb,aupb)
- else
- binorm=anorm_distr(x2,sigma2,alowb,aupb)
- endif
-! print '(a)','Exiting BINORM.'
- return
- end function binorm
-!-----------------------------------------------------------------------------
- real(kind=8) function anorm_distr(x,sigma,alowb,aupb)
-! implicit real*8 (a-h,o-z)
-! to make a normally distributed deviate with zero mean and unit variance
-!
- integer :: iset
- real(kind=8) :: fac,gset,rsq,v1,v2,ran1
- real(kind=8) :: x,sigma,alowb,aupb,gaussdev
- save iset,gset
- data iset/0/
-!elwrite(iout,*) "anorm distr start",x,sigma,alowb,aupb
- if(iset.eq.0) then
- 1 v1=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
- v2=2.0d0*ran_number(0.0d0,1.0d0)-1.0d0
- rsq=v1**2+v2**2
- if(rsq.ge.1.d0.or.rsq.eq.0.0d0) goto 1
- fac=sqrt(-2.0d0*log(rsq)/rsq)
- gset=v1*fac
- gaussdev=v2*fac
- iset=1
- else
- gaussdev=gset
- iset=0
- endif
- anorm_distr=x+gaussdev*sigma
-!elwrite(iout,*) "anorm distr end",x,sigma,alowb,aupb,anorm_distr
- return
- end function anorm_distr
-!-----------------------------------------------------------------------------
- subroutine mult_norm(lda,n,a,x,fail)
-!
-! Generate the vector X whose elements obey the multiple-normal distribution
-! from exp(-0.5*X'AX). LDA is the leading dimension of the moment matrix A,
-! n is the dimension of the problem. FAIL is set at .TRUE., if the smallest
-! eigenvalue of the matrix A is close to 0.
-!
-! implicit double precision (a-h,o-z)
- integer :: lda,n,i,j
- real(kind=8),dimension(lda,n) :: a
- real(kind=8),dimension(n) :: x
- real(kind=8),dimension(100) :: eig,work
- real(kind=8),dimension(3,3) :: vec
- real(kind=8) :: eig_limit=1.0D-8
- logical :: fail
- real(kind=8) :: sigma,alim,xi
- fail=.false.
-! print '(a)','Enter MULT_NORM.'
-!
-! Find the smallest eigenvalue of the matrix A.
-!
-! do i=1,n
-! print '(8f10.5)',(a(i,j),j=1,n)
-! enddo
-#ifdef NAG
- call f02faf('V','U',2,a,lda,eig,work,100,ifail)
-#else
- call djacob(2,lda,10000,1.0d-10,a,vec,eig)
-#endif
-! print '(8f10.5)',(eig(i),i=1,n)
-! print '(a)'
-! do i=1,n
-! print '(8f10.5)',(a(i,j),j=1,n)
-! enddo
- if (eig(1).lt.eig_limit) then
- print *,'From Mult_Norm: Eigenvalues of A are too small.'
- fail=.true.
- return
- endif
-!
-! Generate points following the normal distributions along the principal
-! axes of the moment matrix. Store in WORK.
-!
- do i=1,n
- sigma=1.0D0/dsqrt(eig(i))
- alim=-3.0D0*sigma
- work(i)=anorm_distr(0.0D0,sigma,-alim,alim)
- enddo
-!
-! Transform the vector of normal variables back to the original basis.
-!
- do i=1,n
- xi=0.0D0
- do j=1,n
- xi=xi+a(i,j)*work(j)
- enddo
- x(i)=xi
- enddo
- return
- end subroutine mult_norm
-!-----------------------------------------------------------------------------
- subroutine mult_norm1(lda,n,a,z,box,x,fail)
-!
-! Generate the vector X whose elements obey the multi-gaussian multi-dimensional
-! distribution from sum_{i=1}^m W(i)exp[-0.5*X'(i)A(i)X(i)]. LDA is the
-! leading dimension of the moment matrix A, n is the dimension of the
-! distribution, nlob is the number of lobes. FAIL is set at .TRUE., if the
-! smallest eigenvalue of the matrix A is close to 0.
-!
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
-#ifdef MP
- use MPI_data !include 'COMMON.SETUP'
-#endif
-#ifdef MPI
- include 'mpif.h'
-#endif
- integer :: lda,n,i,j,istep
- real(kind=8),dimension(lda,n) :: a
- real(kind=8),dimension(n) :: z,x
- real(kind=8),dimension(n,n) :: box
- real(kind=8) :: etmp,ww,xi,dec
-! include 'COMMON.IOUNITS'
- logical :: fail
-!
-! Generate points following the normal distributions along the principal
-! axes of the moment matrix. Store in WORK.
-!
-!d print *,'CG Processor',me,' entered MultNorm1.'
-!d print '(2(1pe14.4),3x,1pe14.4)',((a(i,j),j=1,2),z(i),i=1,2)
-!d do i=1,n
-!d print *,i,box(1,i),box(2,i)
-!d enddo
- istep = 0
- 10 istep = istep + 1
- if (istep.gt.10000) then
-! write (iout,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
-! & ' in MultNorm1.'
-! write (*,'(a,i4,2a)') 'CG Processor: ',me,': too many steps',
-! & ' in MultNorm1.'
-! write (iout,*) 'box',box
-! write (iout,*) 'a',a
-! write (iout,*) 'z',z
- fail=.true.
- return
- endif
- do i=1,n
- x(i)=ran_number(box(1,i),box(2,i))
- enddo
- ww=0.0D0
- do i=1,n
- xi=pinorm(x(i)-z(i))
- ww=ww+0.5D0*a(i,i)*xi*xi
- do j=i+1,n
- ww=ww+a(i,j)*xi*pinorm(x(j)-z(j))
- enddo
- enddo
- dec=ran_number(0.0D0,1.0D0)
-! print *,(x(i),i=1,n),ww,dexp(-ww),dec
-!rc if (dec.gt.dexp(-ww)) goto 10
- if(-ww.lt.100) then
- etmp=dexp(-ww)
- else
- return
- endif
- if (dec.gt.etmp) goto 10
-!d print *,'CG Processor',me,' exitting MultNorm1.'
- return
- end subroutine mult_norm1
-!-----------------------------------------------------------------------------
-! ran.f
-!-----------------------------------------------------------------------------
- function ran0(idum)
- integer :: idum
- integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
- integer,parameter :: MASK=123459876
- real(kind=4) :: ran0
- real(kind=4),parameter :: AM=1./IM
- integer :: k
- idum=ieor(idum,MASK)
- k=idum/IQ
- idum=IA*(idum-k*IQ)-IR*k
- if (idum.lt.0) idum=idum+IM
- ran0=AM*idum
- idum=ieor(idum,MASK)
- return
- end function ran0
-! (C) Copr. 1986-92 Numerical Recipes Software *11915
-!-----------------------------------------------------------------------------
- function ran1(idum)
- integer :: idum
- integer,parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
- integer,parameter :: NTAB=32,NDIV=1+(IM-1)/NTAB
- real(kind=4) :: ran1
- real(kind=4),parameter :: AM=1./IM,EPS=1.2e-7,RNMX=1.-EPS
- integer :: j,k,iy
- integer,dimension(NTAB) :: iv
- SAVE iv,iy
- DATA iv /NTAB*0/, iy /0/
- if (idum.le.0.or.iy.eq.0) then
- idum=max(-idum,1)
- do 11 j=NTAB+8,1,-1
- k=idum/IQ
- idum=IA*(idum-k*IQ)-IR*k
- if (idum.lt.0) idum=idum+IM
- if (j.le.NTAB) iv(j)=idum
-11 continue
- iy=iv(1)
- endif
- k=idum/IQ
- idum=IA*(idum-k*IQ)-IR*k
- if (idum.lt.0) idum=idum+IM
- j=1+iy/NDIV
- iy=iv(j)
- iv(j)=idum
- ran1=min(AM*iy,RNMX)
- return
- end function ran1
-! (C) Copr. 1986-92 Numerical Recipes Software *11915
-!-----------------------------------------------------------------------------
- function ran2(idum)
- integer :: idum
- integer,parameter :: IM1=2147483563,IM2=2147483399,IMM1=IM1-1,&
- IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
- NTAB=32,NDIV=1+IMM1/NTAB
- real(kind=4) :: ran2
- real(kind=4),parameter :: AM=1./IM1,EPS=1.2e-7,RNMX=1.-EPS
- integer :: idum2,j,k,iy
- integer,dimension(NTAB) :: iv
- SAVE iv,iy,idum2
- DATA idum2/123456789/, iv/NTAB*0/, iy/0/
- if (idum.le.0) then
- idum=max(-idum,1)
- idum2=idum
- do 11 j=NTAB+8,1,-1
- k=idum/IQ1
- idum=IA1*(idum-k*IQ1)-k*IR1
- if (idum.lt.0) idum=idum+IM1
- if (j.le.NTAB) iv(j)=idum
-11 continue
- iy=iv(1)
- endif
- k=idum/IQ1
- idum=IA1*(idum-k*IQ1)-k*IR1
- if (idum.lt.0) idum=idum+IM1
- k=idum2/IQ2
- idum2=IA2*(idum2-k*IQ2)-k*IR2
- if (idum2.lt.0) idum2=idum2+IM2
- j=1+iy/NDIV
- iy=iv(j)-idum2
- iv(j)=idum
- if(iy.lt.1)iy=iy+IMM1
- ran2=min(AM*iy,RNMX)
- return
- end function ran2
-! (C) Copr. 1986-92 Numerical Recipes Software *11915
-!-----------------------------------------------------------------------------
- function ran3(idum)
- integer :: idum
- integer,parameter :: MBIG=1000000000,MSEED=161803398,MZ=0
-! REAL MBIG,MSEED,MZ
- real(kind=4) :: ran3
- real(kind=4),parameter :: FAC=1./MBIG
-! PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
- integer :: i,iff,ii,inext,inextp,k
- integer :: mj,mk
- integer,dimension(55) :: ma
-! REAL mj,mk,ma(55)
- SAVE iff,inext,inextp,ma
- DATA iff /0/
- if(idum.lt.0.or.iff.eq.0)then
- iff=1
- mj=MSEED-iabs(idum)
- mj=mod(mj,MBIG)
- ma(55)=mj
- mk=1
- do 11 i=1,54
- ii=mod(21*i,55)
- ma(ii)=mk
- mk=mj-mk
- if(mk.lt.MZ)mk=mk+MBIG
- mj=ma(ii)
-11 continue
- do 13 k=1,4
- do 12 i=1,55
- ma(i)=ma(i)-ma(1+mod(i+30,55))
- if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
-12 continue
-13 continue
- inext=0
- inextp=31
- idum=1
- endif
- inext=inext+1
- if(inext.eq.56)inext=1
- inextp=inextp+1
- if(inextp.eq.56)inextp=1
- mj=ma(inext)-ma(inextp)
- if(mj.lt.MZ)mj=mj+MBIG
- ma(inext)=mj
- ran3=mj*FAC
- return
- end function ran3
-! (C) Copr. 1986-92 Numerical Recipes Software *11915
-!-----------------------------------------------------------------------------
-! randgens.f
-!-----------------------------------------------------------------------------
-! $Date: 1994/10/04 16:19:52 $
-! $Revision: 2.1 $
-!
-!
-! See help for RANDOMV on the PSFSHARE disk to understand these
-! subroutines. This is the VS Fortran version of this code.
-!
-!
- subroutine VRND(VEC,N)
-
- use comm_vrandd
-!el integer,dimension(250) :: A
- integer :: LOOP,N !el,I,I147
- integer,dimension(N) :: VEC
-!el COMMON /VRANDD/ A, I, I147
- DO 23000 LOOP=1,N
- I=I+1
- IF(.NOT.(I.GE.251))GOTO 23002
- I=1
-23002 CONTINUE
- I147=I147+1
- IF(.NOT.(I147.GE.251))GOTO 23004
- I147=1
-23004 CONTINUE
- A(I)=IEOR(A(I147),A(I))
- VEC(LOOP)=A(I)
-23000 CONTINUE
- return
- end subroutine VRND
-!-----------------------------------------------------------------------------
- real(kind=8) function RNDV(IDUM)
-
- real(kind=8) :: RM1,RM2,R(99)
- INTEGER :: IA1,IC1,M1, IA2,IC2,M2, IA3,IC3,M3, IDUM
- INTEGER :: IX1,IX2,IX3,J
- SAVE
- DATA IA1,IC1,M1/1279,351762,1664557/
- DATA IA2,IC2,M2/2011,221592,1048583/
- DATA IA3,IC3,M3/15551,6150,29101/
- IF(.NOT.(IDUM.LT.0))GOTO 23006
- IX1 = MOD(-IDUM,M1)
- IX1 = MOD(IA1*IX1+IC1,M1)
- IX2 = MOD(IX1,M2)
- IX1 = MOD(IA1*IX1+IC1,M1)
- IX3 = MOD(IX1,M3)
- RM1 = 1./DBLE(M1)
- RM2 = 1./DBLE(M2)
- DO 23008 J = 1,99
- IX1 = MOD(IA1*IX1+IC1,M1)
- IX2 = MOD(IA2*IX2+IC2,M2)
- R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
-23008 CONTINUE
-23006 CONTINUE
- IX1 = MOD(IA1*IX1+IC1,M1)
- IX2 = MOD(IA2*IX2+IC2,M2)
- IX3 = MOD(IA3*IX3+IC3,M3)
- J = 1+(99*IX3)/M3
- RNDV = R(J)
- R(J) = (DBLE(IX1)+DBLE(IX2)*RM2)*RM1
- IDUM = IX1
- return
- end function RNDV
-!-----------------------------------------------------------------------------
- subroutine VRNDST(SEED)
-
- use comm_vrandd
-!el INTEGER :: A(250)
- INTEGER :: LOOP,IDUM !el,I,I147
- INTEGER :: SEED
-!el real(kind=8) :: RNDV
-!el COMMON /VRANDD/ A, I, I147
- I=0
- I147=103
- IDUM=SEED
- DO 23010 LOOP=1,250
- A(LOOP)=INT(RNDV(IDUM)*2147483647)
-23010 CONTINUE
- return
- end subroutine VRNDST
-!-----------------------------------------------------------------------------
- subroutine VRNDIN(IODEV)
- use comm_vrandd
- INTEGER :: IODEV !el, A(250),I, I147
-!el COMMON/VRANDD/ A, I, I147
- READ(IODEV) A, I, I147
- return
- end subroutine VRNDIN
-!-----------------------------------------------------------------------------
- subroutine VRNDOU(IODEV)
-! This corresponds to VRNDOUT in the APFTN64 version
- use comm_vrandd
- INTEGER :: IODEV !el, A(250),I, I147
-!el COMMON/VRANDD/ A, I, I147
- WRITE(IODEV) A, I, I147
- return
- end subroutine VRNDOU
-!-----------------------------------------------------------------------------
- real(kind=8) function RNUNF(N)
- INTEGER :: IRAN1(2000),N
- real(kind=8) :: FCTOR
- DATA FCTOR /2147483647.0D0/
-! We get only one random number, here! DR 9/1/92
- CALL VRND(IRAN1,1)
- RNUNF= DBLE( IRAN1(1) ) / FCTOR
-!******************************
-! write(6,*) 'rnunf in rnunf = ',rnunf
- return
- end function RNUNF
-!-----------------------------------------------------------------------------
-! readrtns_CSA.F
-!-----------------------------------------------------------------------------
- subroutine random_init(seed)
-!
-! Initialize random number generator
-!
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
-#ifdef MPI
- use MPI_data
- include 'mpif.h'
- logical :: OKRandom !, prng_restart
- real(kind=8) :: r1
- integer :: iseed_array(4)
- integer(kind=8) :: iseed
-#else
- integer :: iseed
-#endif
- real(kind=8) :: seed
- integer :: ierr,error_msg
-! include 'COMMON.IOUNITS'
-! include 'COMMON.TIME1'
-! include 'COMMON.THREAD'
-! include 'COMMON.SBRIDGE'
-! include 'COMMON.CONTROL'
-! include 'COMMON.MCM'
-! include 'COMMON.MAP'
-! include 'COMMON.HEADER'
-! include 'COMMON.CSA'
-! include 'COMMON.CHAIN'
-! include 'COMMON.MUCA'
-! include 'COMMON.MD'
-! include 'COMMON.FFIELD'
-! include 'COMMON.SETUP'
- iseed=-dint(dabs(seed))
- if (iseed.eq.0) then
- write (iout,'(/80(1h*)/20x,a/80(1h*))') &
- 'Random seed undefined. The program will stop.'
- write (*,'(/80(1h*)/20x,a/80(1h*))') &
- 'Random seed undefined. The program will stop.'
-#ifdef MPI
- call mpi_finalize(mpi_comm_world,ierr)
-#endif
- stop 'Bad random seed.'
- endif
-#ifdef MPI
- if (fg_rank.eq.0) then
- seed=seed*(me+1)+1
-#ifdef AMD64
- if(me.eq.king) &
- write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
- OKRandom = prng_restart(me,iseed)
-#else
- do i=1,4
- tmp=65536.0d0**(4-i)
- iseed_array(i) = dint(seed/tmp)
- seed=seed-iseed_array(i)*tmp
- enddo
- if(me.eq.king) &
- write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',&
- (iseed_array(i),i=1,4)
- write (*,*) 'MPI: node= ',me, ' iseed(4)= ',&
- (iseed_array(i),i=1,4)
- OKRandom = prng_restart(me,iseed_array)
-#endif
- if (OKRandom) then
-! r1 = prng_next(me)
- r1=ran_number(0.0D0,1.0D0)
- if(me.eq.king) &
- write (iout,*) 'ran_num',r1
- if (r1.lt.0.0d0) OKRandom=.false.
- endif
- if (.not.OKRandom) then
- write (iout,*) 'PRNG IS NOT WORKING!!!'
- print *,'PRNG IS NOT WORKING!!!'
- if (me.eq.0) then
- call flush(iout)
- call mpi_abort(mpi_comm_world,error_msg,ierr)
- stop
- else
- write (iout,*) 'too many processors for parallel prng'
- write (*,*) 'too many processors for parallel prng'
- call flush(iout)
- stop
- endif
- endif
- endif
-#else
- call vrndst(iseed)
- write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
-#endif
- return
- end subroutine random_init
-!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
- end module random
projects / unres4.git / blobdiff