      SUBROUTINE Haar_USp_2N(n_mezzi,SH,ROUTIN)
C
C     ########################################################################
C
C     Valerio Cappellini                                            April 2007
C
C     USES :: rgnf_lux, GinSE
C     COMPATIBILITY :: Fortran 95
C     VERSION :: 1.0.0
C
C     This program produces 2N x 2N Random Unitary Symplectic Matrices distributed
C     according to their Haar measure, making use of the QR decomposition of 2N x 2N
C     Quaternionic Random (non Self-Dual) Matrices from the Ginibre Ensemble.
C     The QR decomposition is performed by means of (N-1) Householder reflections. 
C
C     Algorithm:
C     ^^^^^^^^^
C      + we start producing a 2N x 2N matrix A(0) = a(0)_{ij} from the Ginibre Ensemble of
C        Quaternionic Random (non Self-Dual) matrices.
C
C                H
C      + we fix S (0) = I_N ( the 2N x 2N identity; here and in the following the 
C        superscript "H" will denote the usual adjoint operation on CPLX matrices)
C
C         + we perform an iterated procedure on the index <k> running from 1 to N-1
C     
C             at each step a matrix H_k is produced and the matrices A(k-1) 
C
C                  H
C             and S (k-1) are upgrated to 
C
C             A(k-1) -> A(k) = H_k * A(k-1) and 
C
C              H          H             H  
C             S (k-1) -> S (k) = H_k * S (k-1)
C
C         + end of iterated procedures
C
C                                     H
C      + We multiply the last raw of S (N-1) by the adjoint of the real quaternion 
C        A(N-1)_{NN} / || A(N-1)_{NN} ||
C
C                    H 
C     Finally SH := S (N-1) is the CPLX Random Unitary Symplectic Matrix distributed 
C     according to the Haar measure and 
C          
C     R = SH * A(0) is an CPLX Upper Triangular Symplectic Matrix with REAL positive 
C     diagonal entries so that this upper positive triangularity could be used as a 
C     test for the subroutine, toghether with the Unitarity of SH.
C     
C     At each step the algorithm make use of the N-k+1 dimensional basis vector 
C
C             / 1  0 | 0  0 | 0  0 | ... ... | 0  0 | 0  0 \ t  , (the superscript t
C     e_k := |       |      |      |         |      |       |      denotes the usual
C             \ 0  1 | 0  0 | 0  0 | ... ... | 0  0 | 0  0 /       transposition) 
C                                    
C     and of the N-k+1 dimensional REAL quaternionic vector
C                                                           t
C     v := (a_{k,k} , a_{k+1,k} , ... , a_{N-1,k} , a_{N,k})   , where each element
C
C                                                      /  z_{ij}   w_{ij}  \
C     a_{ij} is represented by the 2 x 2 CPLX matrix  |   _        _        | .
C                                                      \ -w_{ij}   z_{ij}  /
C
C                                                               H
C     Moreover we use the Euclidean norm || v || = < v | v > = v * v, the first 
C
C                                                                /  z_{kk}   w_{kk}  \
C     entry of the REAL quaternionic vector v , v_1 := a_{kk} = |   _        _        | , 
C                                                                \ -w_{kk}   z_{kk}  /
C     and the phase of the latter P := v_1 / || v_1 || .
C     
C     Then we construct the REAL quaternionic vector u := v + e_k * P * || v || and its 
C
C     positive rescaled-squared-norm c_k := || u ||^2 / 2  = || v || * (|| v || + || v_1 ||)
C     so that finally
C
C
C                /                   |                   \
C               |     I_{2(k-1)}     |         0          |
C               |                    |                    |
C      H_k :=   | ___________________|___________________ |
C               |                    |                    |
C               |                    |                    |
C               |          0         |        M_k         |           
C                \                   |                   /
C     
C     with I_{2(k-1)} being the 2(k-1) x 2(k-1) symplectic identity and M_k the 
C     2(N-k+1) x 2(N-k+1) dimensional CPLX symplectic matrix  
C                                                                 H
C                                     H     /                u * u   \
C     M_k := -(I_{N-k+1} \otimes phase ) * | I_{2(N-k+1)} - --------  |
C                                           \                 c_k    /
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_mezzi,n
      COMPLEX a(2*n_mezzi,2*n_mezzi),SH(2*n_mezzi,2*n_mezzi),
     ,    P(2,2),P_H(2,2)
      REAL ck
      REAL VARIANCE,sum_r
      INTEGER i,j,k

      EXTERNAL ROUTIN

      N=2*n_mezzi
      VARIANCE=1./SQRT(2.)
      call GinSE(n_mezzi,VARIANCE,a,ROUTIN)  ! the matrix A(0) has been produced

      Do i=1,n
          Do j=1,n
              SH(i,j)=(0.,0.)
          Enddo
          SH(i,i)=(1.,0.)
      Enddo                   ! the identity matrix QH(0) has been produced

      sum_r=0.
      do 30 k=1,n-3,2

        do 11 i=k,n
          sum_r=max(sum_r,cabs(a(i,k)))     ! a(k:n,k) contains v
11      continue
        if(sum_r.ne.0.)then                 ! IF v =/= 0
          sum_r=0.
          do 13 i=k,n
            sum_r=sum_r+cabs(a(i,k))**2
13        continue                            ! sum_r = || v ||^2
          P=sqrt(sum_r/                       ! P   := v_1 * || v || / || v_1 ||
     *    (cabs(a(k,k))**2+cabs(a(k,k+1))**2)) 
     *    *a(k:k+1,k:k+1) !                               H
          P_H=CONJG(TRANSPOSE(P))             ! P_H := v_1 * || v || / || v_1 ||
          a(k:k+1,k:k+1)=a(k:k+1,k:k+1)+P     ! a(k:n,k:k+1) now contains u
          P=REAL(MATMUL(a(k:k+1,k:k+1),P_H))  ! P is now diagonal
          ck = (P(1,1)+P(2,2))/2.             ! ck now contains  c_k = ||  u  ||^2 / 2
                                              !  || v || * ( || v || + || v_1 || )
C
C                                                         H
          P_H=-P_H/sqrt(sum_r)                ! P_H = -v_1 / || v_1 ||

          do 16 j=k+2,n-1,2                        
            P = RESHAPE((/(0.,0.)/), (/2, 2/),(/(0.,0.)/))
            do 14 i=k,n-1,2                         
              P=P+MATMUL(CONJG(TRANSPOSE(a(i:i+1,k:k+1))),
     *                   a(i:i+1,j:j+1))
14          continue                            
            P=P/ck                  
            do 15 i=k,n-1,2                     
              a(i:i+1,j:j+1)=a(i:i+1,j:j+1)-MATMUL(a(i:i+1,k:k+1),P)
              a(i:i+1,j:j+1)=MATMUL(P_H,a(i:i+1,j:j+1))
15          continue      
16        continue      
c
c                _N_    _        
c               \   |   u_l * X_{l,j}
c        tau =   |     ---------------  and  X(k:n,j) -> M_k * X(k:n,j)
c               /___|        c_k
c                l=k
c                                   H
c        with X = A (above) or X = S (below)   
c
          do 26 j=1,n-1,2                         
            P = RESHAPE((/(0.,0.)/), (/2, 2/),(/(0.,0.)/))
            do 24 i=k,n-1,2                         
              P=P+MATMUL(CONJG(TRANSPOSE(a(i:i+1,k:k+1))),
     *                   SH(i:i+1,j:j+1))
24          continue                            
            P=P/ck                  
            do 25 i=k,n-1,2                     
              SH(i:i+1,j:j+1)=SH(i:i+1,j:j+1)-MATMUL(a(i:i+1,k:k+1),P)
              SH(i:i+1,j:j+1)=MATMUL(P_H,SH(i:i+1,j:j+1))
25          continue                         
26        continue                          

        endif
30    continue

          P_H=CONJG(TRANSPOSE(a(n-1:n,n-1:n)))/         !  we multiply the 
     /    sqrt(cabs(a(n-1,n-1))**2+cabs(a(n-1,n))**2)   !  last raw of SH(N-1) 
      Do i=1,n-1,2                                      ! 
          SH(n-1:n,i:i+1)= MATMUL(P_H,SH(n-1:n,i:i+1))  !         A(N-1)_{NN} 
      Enddo                                             !  by  -----------------
                                                        !      || A(N-1)_{NN} ||
      return
      END