      SUBROUTINE GUE(N,CMAT,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     A random CPLX Hermitian N x N matrix Z is said to belong to the Gaussian
C     Unitary Ensemble (GUE) if the diagonal elements z_{ii} (which must be real) 
C     and the upper diagonal elements z_{ij} = u_{ij} + i*v_{ij} are independently 
C     chosen with p.d.f.'s
C
C                                   
C                      1         - z_{ii}^2
C     P(z_{ii}) = -----------  e       
C                  SQRT(\pi)
C
C      and
C                                    
C                   2      -2 (u_{ij}^2 + v_{ij}^2)      2      -2 | z_{ij} |^2 
C     P(x_{ij}) = -----  e                           = -----  e
C                  \pi                                  \pi
C
C     respectively. Equivalently
C
C               N(N-1)          N^2
C              --------      - -----
C     P(Z) = 2    2      \pi     2    exp[ - Tr Z* Z ]
C
C     In other words, the diagonal entries have distribution N[0,1/SQRT(2)], where 
C     N[a,b] denotes the "Normal" (or "Gaussian") distribution of mean <a>
C     and variance <b>, while the upper triangular elements have distribution 
C     N[0,1/2] + i*N[0,1/2]
C         
C
C     Z is the output variable <CMAT>.
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
      INTEGER N, emme, mu
      COMPLEX SW(0:1), CMAT(N,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, N
          DO 200 mu = 2*emme+1, 2*N+1
      IF(mu.LE.N+emme+1) CMAT(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 CMAT(emme,mu/2) = CMAT(emme,mu/2) + V / U(1) / 2. * SW(MOD(mu,2))
      IF(mu.eq.2*emme+1) THEN 
          CMAT(emme,mu/2) = CMAT(emme,mu/2) * SQRT(2.)
      ELSEIF(MOD(mu,2).eq.1) THEN
          CMAT(mu/2,emme)=CONJG(CMAT(emme,mu/2))
  200 ENDIF


      RETURN
      END
