module random
!-----------------------------------------------------------------------------
      use io_units
      use prng	! prng.f90 or prng_32.f90
      use math
      implicit none
!
!-----------------------------------------------------------------------------
      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