* questo file contiene alcune routines per la generazione di numeri casuali
* ed un programma test.
* le routines originali si trovano nel file random.f, che contiene anche
* informazioni relative all'autore e riferimenti bibliografici.
* quelle contenute in questo file hanno subito alcune modifiche:
*   1. le variabili sono double precision, anziche' real
*   2. e' possibile chiamare ranmar senza avere preventivamente inizializzato
*      il contatore con rmarin (l'inizializzazione e' automatica)
*   3. e' stata aggiunta la function random() che restituisce un solo valore 
*      casuale double precision
* il programma test illustra come utilizzare le routines per la generazione di
* una matrice di numeri casuali.

*	Giuseppe Rodriguez,  28-04-98
*	Email: rodriguez@unica.it

      program rantest
      integer dim
      parameter (dim=6)
      double precision A(dim,dim)
      integer n, i, j

      n = dim

* rmarin(i1,i2) inizializza una sequenza di numeri casuali.  per cambiare sequenza
* basta modificare i valori di i1 e i2.  se non la si chiama almeno una volta, viene
* inizializzata automaticamente con i1=i2=0.  la seguente chiamata e' quindi ridondante.
      call rmarin(0,0)

* per creare una matrice di numeri casuali la si puo' riempire elemento per elemento
* mediante la funzione random() ...
      do 20 j=1,n
         do 30 i=1,n
	    A(i,j) = random()
   30    continue
   20 continue

      do 10 i=1,n
         print '(1x,10f8.4)', (A(i,j), j=1,n)
   10 continue

      print *

* ... oppure in un colpo solo, chiamando ranmar e passando come parametri l'array da
* riempire ed il numero di elementi dell'array.
      call ranmar(A,n*n)

      do 40 i=1,n
         print '(1x,10f8.4)', (A(i,j), j=1,n)
   40 continue

      end


      subroutine RMARIN(IJ,KL)
C This is the initialization routine for the random number generator RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
C The random number sequences created by these two seeds are of sufficient
C length to complete an entire calculation with. For example, if sveral
C different groups are working on different parts of the same calculation,
C each group could be assigned its own IJ seed. This would leave each group
C with 30000 choices for the second seed. That is to say, this random
C number generator can create 900 million different subsequences -- with
C each subsequence having a length of approximately 10^30.
C
C Use IJ = 1802 & KL = 9373 to test the random number generator. The
C subroutine RANMAR should be used to generate 20000 random numbers.
C Then display the next six random numbers generated multiplied by 4096*4096
C If the random number generator is working properly, the random numbers
C should be:
C           6533892.0  14220222.0  7275067.0
C           6172232.0  8354498.0   10633180.0


      double precision U(97), C, CD, CM
      integer I97, J97
      logical TEST
      common /raset1/ U, C, CD, CM, I97, J97, TEST
      data TEST /.FALSE./

      if( IJ .lt. 0  .or.  IJ .gt. 31328  .or.
     *    KL .lt. 0  .or.  KL .gt. 30081 ) then
          print '(A)', ' The first random number seed must have a value
     *between 0 and 31328'
          print '(A)',' The second seed must have a value between 0 and
     *30081'
            stop
      endif

      i = mod(IJ/177, 177) + 2
      j = mod(IJ    , 177) + 2
      k = mod(KL/169, 178) + 1
      l = mod(KL,     169)

      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
 3       continue
         U(ii) = s
 2    continue

      C = 362436.0 / 16777216.0
      CD = 7654321.0 / 16777216.0
      CM = 16777213.0 /16777216.0

      I97 = 97
      J97 = 33

      TEST = .TRUE.
      return
      end

      subroutine ranmar(RVEC, LEN)
C This is the random number generator proposed by George Marsaglia in
C Florida State University Report: FSU-SCRI-87-50
C It was slightly modified by F. James to produce an array of pseudorandom
C numbers.
      double precision RVEC(*)
      double precision U(97), C, CD, CM
      integer I97, J97
      logical TEST
      common /raset1/ U, C, CD, CM, I97, J97, TEST

      integer ivec

      if( .NOT. TEST ) then
C modified by Giuseppe Rodriguez, rodriguez@unica.it
          call rmarin(0,0)
C         print '(A)',' Call the init routine (RMARIN) before calling RAN
C     *MAR'
C         stop
      endif

      do 100 ivec = 1, LEN
         uni = U(I97) - U(J97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         U(I97) = uni
         I97 = I97 - 1
         if(I97 .eq. 0) I97 = 97
         J97 = J97 - 1
         if(J97 .eq. 0) J97 = 97
         C = C - CD
         if( C .lt. 0.0 ) C = C + CM
         uni = uni - C
         if( uni .lt. 0.0 ) uni = uni + 1.0
         RVEC(ivec) = uni
 100  continue
      return
      end

C added by Giuseppe Rodriguez, rodriguez@unica.it
C this function returns just one double precision random number
      double precision function random()
      double precision x
      call ranmar(x,1)
      random = x
      end
