      SUBROUTINE Haar_O_N(n,Q,ROUTIN)
C
C     ########################################################################
C
C     Valerio Cappellini                                            March 2007
C
C     USES :: rgnf_lux, GinOE
C     COMPATIBILITY :: Fortran 90
C     VERSION :: 1.0.0
C
C     This Program produce N x N REAL Random Orthogonal Matrices distributed according
C     to the Haar measure, making use of the QR decomposition of N x N REAL Random (non 
C     symmetric) Matrices from the Ginibre Ensemble. The QR decomposition is performed
C     by means of (N-1) Householder reflections. 
C
C     Algorithm:
C     ^^^^^^^^^
C      + we start producing an N x N matrix A(0) = a(0)_{ij} from the Ginibre Ensemble.
C      + we fix Q(0) = I_N ( the N x N identity )
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 running matrices A(k-1) 
C             and Q(k-1) are upgrated to A(k-1) -> A(k) = H_k * A(k-1) and 
C             Q(k-1) -> Q(k) = Q(k-1) * H_k
C
C         + end of iterated procedures
C
C      + If A(N-1)_{NN} is negative we change sign to the last column of Q(N-1)
C
C     Finally Q := Q(N-1) is the REAL Random Orthogonal Matrix distributed according 
C     to the Haar measure and 
C
C          t
C     R = Q  * A(0) (here and in the following the superscript t denotes the usual 
C     transposition) is an upper triangular matrix with positive diagonal entries so
C     that this upper positive triangularity could be used as a test for the subroutine,
C     toghether with the orthogonality of Q.
C     
C     At each step the algorithm make use of the N-k+1 dimensional basis vector 
C
C                                    t
C     e_k : (1 , 0 , 0 , ... , 0 , 0)         , of the N-k+1 dimensional vector
C
C                                                           t
C     v := (a_{k,k} , a_{k+1,k} , ... , a_{N-1,k} , a_{N,k})   ,
C
C     its Euclidean norm || v ||  ,  its first entry v_1 = a_{kk}   ,
C
C     the sign of the latter sgn(v_1).
C     
C     Then we construct the vector   u := v + sgn(v_1) * || v || * e_k 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_{k-1}      |         0          |
C               |                    |                    |
C      H_k :=   | ___________________|___________________ |
C               |                    |                    |
C               |                    |                    |
C               |          0         |        M_k         |           
C                \                   |                   /
C     
C     with I_{k-1} being the (k-1) x (k-1) identity , and M_k the (N-k+1) x (N-k+1)
C
C                                                                t
C                                             /             u * u   \
C     dimensional matrix  M_k := -sgn(v_1) * | I_{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
C
      INTEGER n
      REAL a(n,n),Q(n,n),dk,ck
      REAL VARIANCE,sigma,sum,tau
      INTEGER i,j,k

      EXTERNAL ROUTIN

      VARIANCE=1./SQRT(2.)
      call GinOE(n,VARIANCE,a,ROUTIN)  ! the matrix A(0) has been produced

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

      sum=0.
      do 17 k=1,n-1
        do 11 i=k,n
          sum=max(sum,abs(a(i,k)))      ! a(k:n,k) contains v
11      continue
        if(sum.ne.0.)then               ! IF v =/= 0
          sum=0.
          do 13 i=k,n
            sum=sum+a(i,k)**2
13        continue                      ! sum = || v ||^2
          sigma=sign(sqrt(sum),a(k,k))  ! sigma = sgn(v_1) * || v ||
          a(k,k)=a(k,k)+sigma           ! a(k:n,k) now contains u
          ck=sigma*a(k,k)               ! ck now contains c_k = || u ||^2 / 2
                                        !         = || v || * ( || v || + | v_1 | )
          dk=-sign(1.,sigma)            ! dk = - sgn(v_1)

          do 16 j=k+1,n                 !           _N_
            sum=0.                      !          \   |   u_l * a_{l,j}
            do 14 i=k,n                 !   tau =   |     ---------------  
              sum=sum+a(i,k)*a(i,j)     !          /___|        c_k
14          continue                    !           l=k
            tau=sum/ck                  
            do 15 i=k,n                      !
              a(i,j)=dk*(a(i,j)-tau*a(i,k))  !  a(k:n,j) -> M_k * a(k:n,j)
15          continue                         !
16        continue                           

          do 26 i=1,n                   !           _N_
            sum=0.                      !          \   |   Q_{i,l} * u_l 
            do 24 j=k,n                 !   tau =   |     ---------------  
              sum=sum+Q(i,j)*a(j,k)     !          /___|        c_k
24          continue                    !           l=k
            tau=sum/ck                  
            do 25 j=k,n                      !
              Q(i,j)=dk*(Q(i,j)-tau*a(j,k))  !  Q(i,k:n) -> Q(i,k:n) * M_k
25          continue                         !
26        continue

        endif
17    continue
      if(a(n,n).lt.0.) then  ! IF A(N-1)_{NN} is negative THEN
          Do i=1,n           !
              Q(i,n)=-Q(i,n) ! we change sign to the last column of Q(N-1)
          Enddo              !
      Endif                  ! ENDIF
      
      return
      END
