      SUBROUTINE GSE(N,SMAT,ROUTIN)
C
C     ########################################################################
C
C     Valerio Cappellini                                            March 2007
C
C     USES :: rgnf_lux
C     COMPATIBILITY :: Fortran 90
C     VERSION :: 1.0.0
C
C     The Hamiltonian H of a time reversal invariant system must fulfill
C     
C     [T,H] = 0 . The Matrix T, the disrete time reversal symmetry, can either be
C      
C     T^2 = Id     or     T^2 = -Id    and     T=J*K     with       | 0  -1 |
C     K being the operator performing complex conjugation and   J = |       | \otimes I_N
C     The case in which the total angular                           | 1   0 | 
C     momentum is half-odd (odd-spin case, 
C     corresponding to an odd number of spin 1/2 particles) is  T^2 = -1.
C
C     The evolution of physical observables A under Time Reversal (R) is
C
C           R      t  -1                                                R
C     A -> A  = K A  K       (t = transpose) . Obserbles such that A = A  are called
C                                              "self-dual"
C                                      +
C     For Hermitian Hamiltonian ( H = H  ) the self duality condition becomes
C
C                  _  -1   _                                                    +
C            H = K H K    (H = Complex Conjugate of H) that, together with H = H  ,
C
C     means that such self-dual Hermitian matrix H is represented by a Real Quaternion
C     with the following (block) structure:
C
C
C     /  R_{11}     0     |  z_{12}   w_{12}  |  z_{13}   w_{13}  | ... |  z_{1N}   w_{1N}  \  
C    |                    |  _        _       |  _        _       | ... |  _        _        | 
C    |     0      R_{11}  | -w_{12}   z_{12}  | -w_{13}   z_{13}  | ... | -w_{1N}   z_{1N}   |  
C    | ___________________|___________________|___________________| ... |___________________ | 
C    |   _                |                   |                   | ... |                    |  
C    |   z_{12}  -w_{12}  |  R_{22}     0     |  z_{23}   w_{23}  | ... |  z_{2N}   w_{2N}   | 
C    |   _                |                   |  _        _       | ... |  _        _        |  
C    |   w_{12}   z_{12}  |    0      R_{22}  | -w_{23}   z_{23}  | ... | -w_{2N}   z_{2N}   | 
C    | ___________________|___________________|___________________| ... |___________________ |  
C    |   _                |  _                |                   | ... |                    | 
C    |   z_{13}  -w_{13}  |  z_{23}  -w_{23}  |  R_{33}     0     | ... |  z_{3N}   w_{3N}   |  
C    |   _                |  _                |                   | ... |  _        _        | 
C    |   w_{13}   z_{13}  |  w_{23}   z_{23}  |    0      R_{23}  | ... | -w_{3N}   z_{3N}   |  
C    | ___________________|___________________|___________________| ... |___________________ | 
C    |                                                              ...                      |  
C    |  ... ...  ... ...    ... ...  ... ...    ... ...  ... ...  . ... . ... ...  ... ...   | 
C    |  ... ...  ... ...    ... ...  ... ...    ... ...  ... ...  . ... . ... ...  ... ...   |  
C    | ___________________ ___________________ ___________________  ...  ___________________ | 
C    |   _                |  _                |  _                | ... |                    |  
C    |   z_{1N}  -w_{1N}  |  z_{2N}  -w_{2N}  |  z_{3N}  -w_{3N}  | ... |  R_{NN}     0      | 
C    |   _                |  _                |  _                | ... |                    |  
C     \  w_{1N}   z_{1N}  |  w_{2N}   z_{2N}  |  w_{3N}   z_{3N}  | ... |    0      R_{NN}  / 
C 
C
C     The 2N x 2N matrix SMAT(2*N,2*N) contains such a H-matrix.
C
C     We prefere to express H in this CPLX 2N x 2N matrix form... thus each
C     quaternionic element
C                                                   /  z_{ij}   w_{ij}  \
C     q_{ij} is mapped into the 2 x 2 CPLX matrix  |   _        _        |
C                                                   \ -w_{ij}   z_{ij}  /
C
C     A random CPLX Hermitian N x N matrix H with "real quaternionic" elements is 
C     said to belong to the Gaussian Symplectic Ensemble (GSE) if the diagonal elements 
C     q_{ii} of each diagonal real quaternion (which must be real) are independently 
C     chosen with p.d.f.'s
C
C                                   
C                   SQRT(2)      - 2 R_{ii}^2
C     P(R_{ii}) = -----------  e       
C                  SQRT(\pi)
C
C     or equivalently have distribution N[0,1/2], where N[a,b] denotes the "Normal" 
C     (or "Gaussian") distribution of mean <a> and variance <b>.
C
C     The upper triangular off diagonal elements z_{ij} = u_{ij} + i*v_{ij} and
C     w_{ij} = u'_{ij} + i*v'_{ij} are independently chosen with p.d.f.'s
C                                    
C                   4      -4 (u_{ij}^2 + v_{ij}^2)        4      -4 | z_{ij} |^2 
C     P(z_{ij}) = -----  e                             = -----  e
C                  \pi                                    \pi
C
C     and
C                                    
C                   4      -4 (u'_{ij}^2 + v'_{ij}^2)      4      -4 | w_{ij} |^2 
C     P(w_{ij}) = -----  e                             = -----  e
C                  \pi                                    \pi
c
C     or equivalently have distribution N[0,1/2*SQRT(2)] + i*N[0,1/2*SQRT(2)]
C
C     Equivalently
C
C              2*N^2 - 2*N + N/2      -N^2 + N - N/2
C     P(H) = 2                    \pi                 exp[ - Tr H* H ]
C
C
C      For N x N matrices taken from GOE, GUE & GSE 
C     (GSE represantable by 2*N x 2*N CPLX matrices)
C      it holds true:
C
C     =========================+===========================+=============+==============
C                              |                           |             |
C     Kind of elements         | # (indep. REAL variables) |  variance   | Prefactor
C                              |                           |             |
C     =========================+===========================+=============+==============
C                              |                           |             |
C     GOE     diagonal entries |             N             |      1      | 1/SQRT(2\pi)
C                              |                           |             |
C     -------------------------+---------------------------+-------------+--------------
C                              |                           |             |
C     GOE off-diagonal entries |          N(N-1)/2         |  1/SQRT(2)  | 1/SQRT(\pi)
C                              |                           |             |
C     -------------------------+---------------------------+-------------+--------------
C                              |                           |             |
C     GUE     diagonal entries |             N             |  1/SQRT(2)  | 1/SQRT(\pi)
C                              |                           |             |
C     -------------------------+---------------------------+-------------+--------------
C                              |                           |             |
C     GUE off-diagonal entries |           N(N-1)          |     1/2     | SQRT(2/\pi)
C                              |                           |             |
C     -------------------------+---------------------------+-------------+--------------
C                              |                           |             |
C     GSE     diagonal entries |             N             |     1/2     | SQRT(2/\pi)
C                              |                           |             |
C     -------------------------+---------------------------+-------------+--------------
C                              |                           |             |
C     GSE off-diagonal entries |         2*N(N-1)          | 1/2*SQRT(2) | 2/SQRT(\pi)
C                              |                           |             |
C     =========================+===========================+=============+==============
C
C
C     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     ROUTIN   : Name of the subroutine you are going to use for generating
C                the random numbers (suggested :: rgnf_lux.f)
C
C     IMPORTANT [1] !!! Remember to initialize the Subroutine rgnf_lux.f
C                       by means of "Call Ranset_()"
C     IMPORTANT [2] !!! Remember to put the file "i_seed_rng" into the 
C                       working directory.
C     IMPORTANT [3] !!! Remember to insert main_program: "EXTERNAL rgnf_lux"
C
C     ########################################################################
C
$     debug
C
      INTEGER N, emme, mu
      COMPLEX SW(0:1), SMAT(2*N,2*N)
      REAL U(2), S, T, A, B, R1, R2, V, X, Y, Q
      SAVE  S, T, A, B, R1, R2, SW
      EXTERNAL ROUTIN
      DATA  S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/
      DATA  R1, R2 ,SW / 0.27597, 0.27846, (0. , 1.), (1. , 0.)/
C         generate pair of uniform deviates

      DO 200 emme = 1, 2*N-1,2
          DO 200 mu = 2*emme+1, 4*N+1
      IF(mu.LE.2*N+emme+1) SMAT(emme,mu-emme-1) = (0.,0.)
   50 CALL ROUTIN(U,2)
      V = 1.7156 * (U(2) - 0.5)
      X = U(1) - S
      Y = ABS(V) - T
      Q = X**2 + Y*(A*Y - B*X)
C           accept P if inside inner ellipse
      IF (Q .LT. R1)  GO TO 100
C           reject P if outside outer ellipse
      IF (Q .GT. R2)  GO TO 50
C           reject P if outside acceptance region
      IF (V**2 .GT. -4.0 *ALOG(U(1)) *U(1)**2)  GO TO 50
C           ratio of P's coordinates is normal deviate
  100 SMAT(emme,mu/2)=SMAT(emme,mu/2)+V/U(1)/2./SQRT(2.)*SW(MOD(mu,2))
      IF(mu.eq.2*emme+1) THEN
          SMAT(emme  ,mu/2  ) =        SMAT(emme,mu/2) * SQRT(2.)
          SMAT(emme+1,mu/2+1) =  CONJG(SMAT(emme,mu/2))
      ELSEIF(mu.eq.2*emme+3) THEN
          SMAT(emme  ,mu/2  ) = (0.,0.)
          SMAT(emme+1,mu/2-1) = (0.,0.)
      ELSEIF(MOD(mu,4).eq.1) THEN
          SMAT(emme+1,mu/2  ) =  CONJG(SMAT(emme,mu/2-1))
          SMAT(emme+1,mu/2-1) = -CONJG(SMAT(emme,mu/2  ))
      ENDIF

  200 Continue

      DO 300 emme = 1, 2*N
          DO 300 mu = 1, emme-1
          SMAT(emme,mu) =  CONJG(SMAT(mu,emme))
  300 Continue

      RETURN
      END