! Main program for structured eigenvalue computation.
!
! This program finds the roots of a given polynomial
! by computing the eigenvalues of the associated
! companion pencil with O(n^2) complexity. 
! The unitary + rank one and quasiseparable
! structures of the pencil are used.
!
! This is an implementation
! of the method described in the preprint
! "Implicit QR for Companion-Like Pencils",
! by P. Boito, Y. Eidelman and L. Gemignani.
!
! Single shift only for now.
!
! Lapack needs to be installed in order to evaluate
! accuracy of computed eigenvalues.
! The Fortran routine ssort.f is also needed 
! (it can easily be found on the Internet). 
!
! Last update on Sept 17, 2014.
!
! Please send comments or bug reports to
! paola.boito@unilim.fr
!
! (Edit: slight modification to givens subroutine on Dec 19, 2016).


program main

implicit none

integer :: n,i,j,k,lwork,info,lapacksize 
integer :: clock_start,clock_end,clock_rate
real(kind(1d0)) :: w1,w2,w3,our_time,maxdiffeR,maxdiffeI,lapack_time,normpoly
real(kind(1d0)), allocatable :: rwork(:),dist1(:),dist2(:) 
real(kind(1d0)), allocatable :: eR1(:),eR2(:),eI1(:),eI2(:),seR1(:),seR2(:),seI1(:),seI2(:),diffeR(:),diffeI(:)
complex(kind(1d0)), allocatable :: work(:),vl(:,:),vr(:,:) ,alphaL(:),betaL(:)
complex(kind(1d0)), allocatable :: eig(:),myeig(:),Am(:,:),Bm(:,:),poly(:)

write(*,"(a)",advance="no") "Degree = " 
read(*,*) n 

allocate(poly(n+1),myeig(n)) 

! Define polynomial coefficients.
! (Comment/uncomment as needed).
!
!
!*******************************************  
! TEST 1: Random coefficients 
!*******************************************
   
call random_seed
   
do i=1,n+1 
   call random_number(w1) 
   call random_number(w2) 
   poly(i)=(2.d0,0.d0)*w1-(1.d0, 0.d0)+(0.d0,1.d0)*w2-(0.d0, 1.d0) 
enddo 

!
!********************************************
! TEST 2a: Cyclotomic polynomial  (roots of 1)
!********************************************
!
!poly(1)=1.d0
!poly(2:n)=0.d0
!poly(n+1)=(-1.d0,0.d0)

!********************************************
! TEST 2b: Cyclotomic polynomial (roots of i)
!********************************************

!poly(1)=1.d0
!poly(2:n)=0.d0
!poly(n+1)=(0.d0,-1.d0)

!********************************************  
! TEST 3: From file 
!********************************************
!
!open(unit=300,file="mypoly.txt")
!do i=1,n+1
!  write(300,*) poly(i)
!end do
!close(300)
!
!****************************************  
! TEST 4: All coefficients = 1
!****************************************
!
!do i=1,n+1
!   poly(i)=1.d0
!enddo
!
!****************************************

! Normalize polynomial:
normpoly=0.d0
do i=1,n+1
   normpoly=normpoly+poly(i)*conjg(poly(i))
enddo
normpoly=sqrt(normpoly)
poly=poly/normpoly   


! Apply fast_companion_QZ.f90 to poly and compute running time


  CALL SYSTEM_CLOCK(COUNT_RATE=clock_rate) 
  CALL SYSTEM_CLOCK(COUNT=clock_start) 

  call fast_companion_QZ(n,poly,myeig)

  CALL SYSTEM_CLOCK(COUNT=clock_end)  
  our_time=(REAL(clock_end - clock_start)/clock_rate) 
  write(*,*)"Our time = ",our_time 


! Apply the Lapack routine ZGEGV to the companion pencil (Am,Bm)
! associated with POLY and compute running time

allocate(vL(1,n),vR(1,n),Am(n,n),Bm(n,n),eig(n),dist1(n),dist2(n),alphaL(n),betaL(n))    

Am=0.d0 
do i=1,n-1 
   Am(i+1,i)=(1.d0,0.d0) 
enddo 
do i=1,n
   Am(1,i)=-poly(i+1) 
enddo   
Bm=0.d0
Bm(1,1)=poly(1)
do i=2,n
   Bm(i,i)=1.d0
enddo   


  LWORK = max(1,6*n+100) 
  allocate (WORK(LWORK)) 
  allocate(rwork(8*n)) 
  CALL SYSTEM_CLOCK(COUNT=clock_start)
  call ZGEGV('n','n',n,Am,n,Bm,n,alphaL,betaL,vL,1,vR,1,WORK,LWORK,RWORK,INFO) 

  CALL SYSTEM_CLOCK(COUNT=clock_end)  
  lapack_time=(REAL(clock_end - clock_start)/clock_rate) 
  write(*,*)"Lapack time = ",lapack_time


! Compare eigenvalues:

do j=1,n
   eig(j)=alphaL(j)/betaL(j)
enddo
  !
  ! comparison by minimum distance  

  compare_outer:do j=1,n 
     compare_inner:do i=1,n 
        dist2(i)=abs(myeig(j)-eig(i)) 
     end do compare_inner 
     dist1(j)=minval(dist2) 
  end do compare_outer 
  write(*,*)"distance = ",maxval(dist1) 
   

  ! comparison by sorting
  ! (warning: this might sometimes be unreliable, due to complex conjugate
  ! eigenvalues)
  allocate(eR1(n),eR2(n),eI1(n),eI2(n),seR1(n),seR2(n),seI1(n),seI2(n),diffeR(n),diffeI(n))

  eR1=real(myeig)
  eI1=aimag(myeig)
  eR2=real(eig)
  eI2=aimag(eig)

  call sort_eig(n,eR1,eI1,seR1,seI1)
  call sort_eig(n,eR2,eI2,seR2,seI2)

  diffeR = abs(seR2 - seR1)
  diffeI = abs(seI2 - seI1)
  maxdiffeR = maxval(abs(diffeR))
  maxdiffeI = maxval(abs(diffeI))

  write(*,*)"distance (by sorting) = ",maxdiffeR, maxdiffeI
  !write(*,*)"abs. distance (by sorting) = ",sqrt(maxdiffeR**2 + maxdiffeI**2)

!write(*,*)"myeig=",myeig
!write(*,*)"Lapack eig=",eig
!write(*,*)"Lapack alpha=",alphaL
!write(*,*)"Lapack beta=",betaL

end program main



!***********************************************************


subroutine fast_companion_QZ(n,poly,myeig)

! Fast eigenvalue computation.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: poly(n+1)
complex(kind(1d0)), intent(out) :: myeig(n)


integer :: i,j,k,kin,kfin,numit,total_it,Bdefl
real(kind(1d0)) :: tol,eps,ttt,w1rand,w2rand
complex(kind(1d0)) :: adiag(n),bdiag(n)
complex(kind(1d0)) :: sigmaa(n-1),gv(n,2),hv(2,n),bv(2,2,n-1),z(n),w(n)
complex(kind(1d0)) :: db(n),gu(n-1),hu(n),bu(n-1),p(n),q(n)
complex(kind(1d0)) :: xx(2),G(2,2),tkin,tkfin,tj,tjm1,yy(2)

! Variables from structured QZ subroutine:

complex(kind(1d0)) :: sigmaanew(n-1),gvnew(n,2),hvnew(2,n),bvnew(2,2,n-1)
complex(kind(1d0)) :: gunew(n-1),hunew(n),bunew(n-1),Qin(2,2)
complex(kind(1d0)) :: du(n),du1tilde(n+1),gu1tilde(n+1,2),hu1tilde(2,n+1),bu1tilde(2,2,n+1),w1(n),z1(n)
complex(kind(1d0)) :: dv1tilde(n+1),gv1tilde(n+1,3),hv1tilde(3,n+2),bv1tilde(3,3,n+1),sigmav(n-1),p1(n),q1(n)
complex(kind(1d0)) :: db1(n),du1(n),gu1(n,2),hu1(2,n),bu1(2,2,n),dv1(n),gv1(n,3),hv1(3,n),bv1(3,3,n)

complex(kind(1d0)) :: QQ(2,2),epsi,auxmat(2,3),GGamma(2,4),fa,fv,phia(2),phiv(2),auxvec(2),chi
complex(kind(1d0)) :: gggamma,c,theta,fu,phiu,fb,phi,vareps,PPhi(2,2),auxmat2(2,2),Omega(2,2),ZZ(2,2)
complex(kind(1d0)) :: newQ(2,2),auxmat24(2,4),auxmat32(3,2),Gammarow(2),auxmat22(2,2),htilde1,betav(3)
complex(kind(1d0)) :: bs(2,2,n-1),gs(n,2),hs(2,n),Sold(2,3),Snew(2,3),yold,ynew,zvec(2),VV(2,2),AA(3,2)
complex(kind(1d0)) :: auxvec3(3),FF(3,3),WW(2,2),xnum,yvec(2),znum,Effe(2,2),Zmat(2,2),CC(3,4),Gammap(2,4)
complex(kind(1d0)) :: Lambda(2,3),betau(2),hutilde,hvtilde1,GG(2,2),auxmat34(3,4)

! *** for checks only *************************************************
! complex(kind(1d0)) :: A(n,n),B(n,n),V(n,n),U(n,n),mat2n(2,n),matn2(n,2),othera(n,n),otherb(n,n),normaF2,normbF2
! *********************************************************************


eps=epsilon(1d0)
tol=10*eps
total_it=0


! Compute initial quasiseparable generators

call companion_generators(n,poly,sigmaa,gv,hv,bv,z,w,db,gu,hu,bu,p,q)

! ******* CHECK ******************
!call recover_A(A,sigmaa,gv,hv,bv,z,w,n)
!call recover_B(B,db,gu,hu,bu,p,q,n)
!V=0.d0
!do j=1,n-1
!   V(j+1,j)=1.d0   
!enddo
!V(1,n)=1.d0
!
!U=0.d0
!do j=1,n
!   U(j,j)=1.d0
!enddo
!
! *********************************


! Perform deflation on A and on B:
! determine the indices kin and kfin 
! where the QZ process should start and then stop

kfin=n
numit=0

do while (kfin>1)
   Bdefl=1
   do while (Bdefl==1)
      k=kfin-1
      do while (k>0 .AND. abs(db(k))>tol)
         k=k-1
      enddo   
      if (k>0) then
         write(*,*) "Deflation on B is needed."
         return
      else
      Bdefl=0   
      end if
   enddo
   ! Deflation on A: choose kfin
   j=kfin   
   tj=gv(j,1)*hv(1,j)+gv(j,2)*hv(2,j)
   tjm1=gv(j-1,1)*hv(1,j-1)+gv(j-1,2)*hv(2,j-1)
   do while (j>1 .AND. abs(sigmaa(j-1))<eps*(abs(tj-z(j)*conjg(w(j)))+abs(tjm1-z(j-1)*conjg(w(j-1)))))   
      !write(*,*) "deflation"
      adiag(j)=gv(j,1)*hv(1,j)+gv(j,2)*hv(2,j)-z(j)*conjg(w(j))
      bdiag(j)=db(j)
      numit=0
      j=j-1
   enddo
   kfin=j
   ! Choose kin with single subdiag criterion
   kin=kfin-1
   tkin=gv(kin,1)*hv(1,kin)+gv(kin,2)*hv(2,kin)
   ttt=abs(tkin-z(kin)*conjg(w(kin)))
   do while (kin>1 .AND. abs(sigmaa(kin-1))>eps*ttt)
      kin=kin-1
   enddo
   ! Now run the QZ algorithm
   if (kin>0) then
      ! Do the "magic trick" on A, increase kin if necessary
      j=kfin-1
      if (numit<10) then
         tj=gv(j,1)*hv(1,j)+gv(j,2)*hv(2,j)
         tkfin=gv(kfin,1)*hv(1,kfin)+gv(kfin,2)*hv(2,kfin)
         xx(1)=db(kfin)*(tj-z(j)*conjg(w(j)))-(tkfin-z(kfin)*conjg(w(kfin)))*db(j)
         xx(2)=db(kfin)*sigmaa(j)
      else
         !write(*,*) "Exceptional shift"   
         call random_number(w1rand) 
         call random_number(w2rand) 
         xx(1)=(2.d0,0.d0)*w1rand-(1.d0, 0.d0)+(0.d0,1.d0)*w2rand-(0.d0, 1.d0) 
         call random_number(w1rand) 
         call random_number(w2rand) 
         xx(2)=(2.d0,0.d0)*w1rand-(1.d0, 0.d0)+(0.d0,1.d0)*w2rand-(0.d0, 1.d0) 
      end if   
      call givens(xx,G)
      if (j>1) then ! do the "magic trick", otherwise we already know kin=1
         yy(1)=sigmaa(j-1)*G(1,1)
         yy(2)=sigmaa(j-1)*G(2,1)
         do while (j>kin .AND. abs(yy(2))>tol)
            j=j-1
            if (numit<10) then
               ! Choose shift
               tj=gv(j,1)*hv(1,j)+gv(j,2)*hv(2,j)
               tkfin=gv(kfin,1)*hv(1,kfin)+gv(kfin,2)*hv(2,kfin)
               xx(1)=db(kfin)*(tj-z(j)*conjg(w(j)))-(tkfin-z(kfin)*conjg(w(kfin)))*db(j)
               xx(2)=db(kfin)*sigmaa(j)
               call givens(xx,G)            
            end if
            if (j>1) then
               yy=sigmaa(j-1)*G(:,1)
            end if
         enddo
         
         ! Update kin
         kin=j
      else      
         kin=1
      end if
      ! Reinitialize counter after exceptional shift
         if (numit>=10) then
            numit=0
         end if

      
      
! Perform structured QZ from kin to kfin  
   
!*************************************
!*************************************    
     

! Recall that:
! - sigma,gv,hv,bv,z,w are parameters for the matrix A,
! - d,gu,hu,bu,p,q are parameters for the matrix B,
! - Qin is the initial Givens matrix (encoding the shift).

Qin=G

! In order to speed up the program:
! if kfin is large do the loop until the end,
! otherwise stop early
!*******************************************
if (kfin>n-6) then
!*******************************************

do j=1,n-1
   sigmav(j)=sigmaa(j)+z(j+1)*conjg(w(j))
enddo   
epsi=gv(1,1)*hv(1,1)+gv(1,2)*hv(2,1)-z(1)*conjg(w(1))
if (kin==1) then
    QQ=Qin
else
    QQ(1,1)=1.d0
    QQ(2,2)=1.d0
    QQ(1,2)=0.d0
    QQ(2,1)=0.d0
endif

! *** (1) ****

GGamma=0.d0
auxmat(1,1)=gv(1,1)*hv(1,1)+gv(1,2)*hv(2,1)
auxmat(1,2)=gv(1,1)*bv(1,1,1)+gv(1,2)*bv(2,1,1)
auxmat(1,3)=gv(1,1)*bv(1,2,1)+gv(1,2)*bv(2,2,1)
auxmat(2,1)=sigmav(1)
auxmat(2,2)=gv(2,1)
auxmat(2,3)=gv(2,2)
GGamma(:,1:3)=matmul(QQ,auxmat)

gv1tilde(2,:)=GGamma(1,1:3)
fv=GGamma(2,1)
phiv=GGamma(2,2:3)
auxvec(1)=QQ(1,1)*z(1)+QQ(1,2)*z(2)
auxvec(2)=QQ(2,1)*z(1)+QQ(2,2)*z(2)

z1(1)=auxvec(1)
chi=auxvec(2)
fa=fv-chi*conjg(w(1))
phia=phiv

! ******* CHECK ******************
!mat2n=A(1:2,:)
!A(1:2,:)=matmul(QQ,mat2n)
!mat2n=V(1:2,:)
!V(1:2,:)=matmul(QQ,mat2n)
! ********************************


! *** (2) ****

gggamma=w(1)
c=p(1)
theta=q(1)
du(1)=db(1)+p(1)*conjg(q(1))
fu=du(1)
phiu=gu(1)
fb=db(1)
phi=gu(1)

! *** (3) ****

do k=1,n-2
   ! *** (3a) ****
   vareps=phi*hu(k+1)-c*conjg(q(k+1))
   epsi=phia(1)*hv(1,k+1)+phia(2)*hv(2,k+1)-chi*conjg(w(k+1))
   du(k+1)=db(k+1)+p(k+1)*conjg(q(k+1))
   auxmat2(1,1)=fb
   auxmat2(1,2)=vareps
   auxmat2(2,1)=0.d0
   auxmat2(2,2)=db(k+1)
   PPhi=matmul(QQ,auxmat2)
   
   
   ! ******* CHECK ******************
   !mat2n=B(k:k+1,:)
   !B(k:k+1,:)=matmul(QQ,mat2n)
   !mat2n=U(k:k+1,:)
   !U(k:k+1,:)=matmul(QQ,mat2n)
   ! ********************************
   
   
   ! *** (3b) ****
   if (k>=kin) then
      xx(1)=conjg(PPhi(2,1))
      xx(2)=conjg(PPhi(2,2))
      call givens(xx,GG)
      ZZ(1,1)=conjg(GG(2,1))
      ZZ(1,2)=conjg(GG(1,1))
      ZZ(2,1)=conjg(GG(2,2))
      ZZ(2,2)=conjg(GG(1,2))
   else
      ZZ(1,1)=1.d0
      ZZ(1,2)=0.d0
      ZZ(2,1)=0.d0
      ZZ(2,2)=1.d0
   endif
   auxmat2(1,1)=fa
   auxmat2(1,2)=epsi
   auxmat2(2,1)=0.d0
   auxmat2(2,2)=sigmaa(k+1)
   
   Omega=matmul(auxmat2,ZZ)   
   
   ! ****** CHECK ********
   !matn2=A(:,k:k+1)
   !A(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=V(:,k:k+1)
   !V(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=B(:,k:k+1)
   !B(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=U(:,k:k+1)
   !U(:,k:k+1)=matmul(matn2,ZZ)
   ! **********************
   
   
   ! *** (3c) ****
   if (k>=kin .AND. k<=kfin-2) then
      xx=Omega(:,1)
      call givens(xx,newQ)
      sigmaanew(k)=newQ(1,1)*xx(1)+newQ(1,2)*xx(2)
   else
      if (k+1==kin) then
         newQ=Qin
      else
         newQ(1,1)=1.d0
         newQ(1,2)=0.d0
         newQ(2,1)=0.d0
         newQ(2,2)=1.d0                     
      end if
      sigmaanew(k)=newQ(1,1)*Omega(1,1)+newQ(1,2)*Omega(2,1)
   endif
   
   
   ! ***** CHECK *********
   !mat2n=A(k+1:k+2,:)
   !A(k+1:k+2,:)=matmul(newQ,mat2n)
   !mat2n=V(k+1:k+2,:)
   !V(k+1:k+2,:)=matmul(newQ,mat2n)
   ! *********************
   
   
   
   ! *** (3d) ****
   auxmat24(1,1)=fv
   auxmat24(1,2)=phiv(1)*hv(1,k+1)+phiv(2)*hv(2,k+1)
   auxmat24(1,3)=phiv(1)*bv(1,1,k+1)+phiv(2)*bv(2,1,k+1)
   auxmat24(1,4)=phiv(1)*bv(1,2,k+1)+phiv(2)*bv(2,2,k+1)
   auxmat24(2,1)=z(k+2)*conjg(gggamma)
   auxmat24(2,2)=sigmav(k+1)
   auxmat24(2,3)=gv(k+2,1)
   auxmat24(2,4)=gv(k+2,2)
   Gammap=matmul(newQ,auxmat24)
   GGamma=Gammap
   auxmat2(1,1)=Gammap(1,1)*ZZ(1,1)+Gammap(1,2)*ZZ(2,1)
   auxmat2(1,2)=Gammap(1,1)*ZZ(1,2)+Gammap(1,2)*ZZ(2,2)
   auxmat2(2,1)=Gammap(2,1)*ZZ(1,1)+Gammap(2,2)*ZZ(2,1)
   auxmat2(2,2)=Gammap(2,1)*ZZ(1,2)+Gammap(2,2)*ZZ(2,2)
   
   GGamma(1,1)=auxmat2(1,1)
   GGamma(1,2)=auxmat2(1,2)
   GGamma(2,1)=auxmat2(2,1)
   GGamma(2,2)=auxmat2(2,2)
   auxmat34(1,1)=1.d0
   auxmat34(1,2)=0.d0
   auxmat34(1,3)=0.d0
   auxmat34(1,4)=0.d0
   auxmat34(2,1)=0.d0
   auxmat34(3,1)=0.d0
   auxmat34(2:3,2)=hv(:,k+1)
   auxmat34(2:3,3:4)=bv(:,:,k+1)
   CC=auxmat34
   auxmat32=matmul(auxmat34(:,1:2),ZZ)
   CC(1,1)=auxmat32(1,1)
   CC(1,2)=auxmat32(1,2)
   CC(2,1)=auxmat32(2,1)
   CC(2,2)=auxmat32(2,2)
   CC(3,1)=auxmat32(3,1)
   CC(3,2)=auxmat32(3,2)
   dv1tilde(k+2)=GGamma(1,1)
   gv1tilde(k+2,:)=GGamma(1,2:4)
   betav=GGamma(2,2:4)
   fv=betav(1)
   phiv=betav(2:3)
   hv1tilde(1,k+2)=CC(1,1)
   hv1tilde(2,k+2)=CC(2,1)
   hv1tilde(3,k+2)=CC(3,1)
   bv1tilde(1,1,k+2)=CC(1,2)
   bv1tilde(1,2,k+2)=CC(1,3)
   bv1tilde(1,3,k+2)=CC(1,4)
   bv1tilde(2,1,k+2)=CC(2,2)
   bv1tilde(2,2,k+2)=CC(2,3)
   bv1tilde(2,3,k+2)=CC(2,4)
   bv1tilde(3,1,k+2)=CC(3,2)
   bv1tilde(3,2,k+2)=CC(3,3)
   bv1tilde(3,3,k+2)=CC(3,4)
   ! *** (3e) ****
   auxmat(1,1)=fu
   auxmat(1,2)=phiu*hu(k+1)
   auxmat(1,3)=phiu*bu(k+1)
   auxmat(2,1)=p(k+1)*conjg(theta)
   auxmat(2,2)=du(k+1)
   auxmat(2,3)=gu(k+1)
   auxmat2=matmul(auxmat(:,1:2),ZZ)
   auxmat(1,1)=auxmat2(1,1)
   auxmat(1,2)=auxmat2(1,2)
   auxmat(2,1)=auxmat2(2,1)
   auxmat(2,2)=auxmat2(2,2)
   Lambda=matmul(QQ,auxmat)
   auxmat(1,1)=1.d0
   auxmat(1,2)=0.d0
   auxmat(1,3)=0.d0
   auxmat(2,1)=0.d0
   auxmat(2,2)=hu(k+1)
   auxmat(2,3)=bu(k+1)
   auxmat2=matmul(auxmat(:,1:2),ZZ)
   auxmat(1,1)=auxmat2(1,1)
   auxmat(1,2)=auxmat2(1,2)
   auxmat(2,1)=auxmat2(2,1)
   auxmat(2,2)=auxmat2(2,2)
   du1tilde(k)=Lambda(1,1)
   gu1tilde(k,1)=Lambda(1,2)
   gu1tilde(k,2)=Lambda(1,3)
   betau=Lambda(2,2:3)
   fu=betau(1)
   phiu=betau(2)
   hu1tilde(1,k+1)=auxmat(1,1)
   hu1tilde(2,k+1)=auxmat(2,1)
   bu1tilde(1,1,k+1)=auxmat(1,2)
   bu1tilde(1,2,k+1)=auxmat(1,3)
   bu1tilde(2,1,k+1)=auxmat(2,2)
   bu1tilde(2,2,k+1)=auxmat(2,3)
   ! *** (3f) ****
   auxvec(1)=newQ(1,1)*chi+newQ(1,2)*z(k+2)
   auxvec(2)=newQ(2,1)*chi+newQ(2,2)*z(k+2)
   z1(k+1)=auxvec(1)
   chi=auxvec(2)
   auxvec(1)=conjg(ZZ(1,1))*gggamma+conjg(ZZ(2,1))*w(k+1)
   auxvec(2)=conjg(ZZ(1,2))*gggamma+conjg(ZZ(2,2))*w(k+1)   
   w1(k)=auxvec(1)
   gggamma=auxvec(2)
   ! *** (3g) ****
   auxvec(1)=conjg(ZZ(1,1))*theta+conjg(ZZ(2,1))*q(k+1)
   auxvec(2)=conjg(ZZ(1,2))*theta+conjg(ZZ(2,2))*q(k+1)
   q1(k)=auxvec(1)
   theta=auxvec(2)
   auxvec(1)=QQ(1,1)*c+QQ(1,2)*p(k+1)
   auxvec(2)=QQ(2,1)*c+QQ(2,2)*p(k+1)
   p1(k)=auxvec(1)
   c=auxvec(2)
   ! *** (3h) ****
   fa=fv-chi*conjg(gggamma)
   phia=phiv
   fb=fu-c*conjg(theta)
   phi=phiu
   QQ=newQ  
   
enddo

   

! *** (4) ****

vareps=phi*hu(n)-c*conjg(q(n))
du(n)=db(n)+p(n)*conjg(q(n))
auxmat2(1,1)=fb
auxmat2(1,2)=vareps
auxmat2(2,1)=0.d0
auxmat2(2,2)=db(n)
PPhi=matmul(QQ,auxmat2)

! ***** CHECK ******
!mat2n=B(n-1:n,:)
!B(n-1:n,:)=matmul(QQ,mat2n)
!mat2n=U(n-1:n,:)
!U(n-1:n,:)=matmul(QQ,mat2n)
! ******************


! *** (5) ****

xx(1)=conjg(PPhi(2,1))
xx(2)=conjg(PPhi(2,2))
call givens(xx,GG)
ZZ(1,1)=conjg(GG(2,1))
ZZ(1,2)=conjg(GG(1,1))
ZZ(2,1)=conjg(GG(2,2))
ZZ(2,2)=conjg(GG(1,2))

! *** (6) ****

auxvec(1)=fv
auxvec(2)=phiv(1)*hv(1,n)+phiv(2)*hv(2,n)
Gammarow(1)=auxvec(1)*ZZ(1,1)+auxvec(2)*ZZ(2,1)
Gammarow(2)=auxvec(1)*ZZ(1,2)+auxvec(2)*ZZ(2,2)
auxmat32(1,1)=1.d0
auxmat32(1,2)=0.d0
auxmat32(2,1)=0.d0
auxmat32(2,2)=hv(1,n)
auxmat32(3,1)=0.d0
auxmat32(3,2)=hv(2,n)
CC(:,1:2)=matmul(auxmat32,ZZ)
hvtilde1=1.d0
dv1tilde(n+1)=Gammarow(1)
gv1tilde(n+1,1)=Gammarow(2)
gv1tilde(n+1,2)=0.d0
gv1tilde(n+1,3)=0.d0
hv1tilde(1,n+1)=CC(1,1)
hv1tilde(2,n+1)=CC(2,1)
hv1tilde(3,n+1)=CC(3,1)
bv1tilde(1,1,n+1)=CC(1,2)
bv1tilde(2,1,n+1)=CC(2,2)
bv1tilde(3,1,n+1)=CC(3,2)
bv1tilde(1,2,n+1)=0.d0
bv1tilde(2,2,n+1)=0.d0
bv1tilde(3,2,n+1)=0.d0
bv1tilde(1,3,n+1)=0.d0
bv1tilde(2,3,n+1)=0.d0
bv1tilde(3,3,n+1)=0.d0


! ***** CHECK *****
!matn2=A(:,n-1:n)
!A(:,n-1:n)=matmul(matn2,ZZ)
!matn2=V(:,n-1:n)
!V(:,n-1:n)=matmul(matn2,ZZ)
!matn2=B(:,n-1:n)
!B(:,n-1:n)=matmul(matn2,ZZ)
!matn2=U(:,n-1:n)
!U(:,n-1:n)=matmul(matn2,ZZ)
! *****************




! *** (7) ****

auxmat2(1,1)=fu
auxmat2(1,2)=phiu*hu(n)
auxmat2(2,1)=p(n)*conjg(theta)
auxmat2(2,2)=du(n)
auxmat22=matmul(auxmat2,ZZ)
Lambda(:,1:2)=matmul(QQ,auxmat22)
auxmat2(1,1)=1.d0
auxmat2(1,2)=0.d0
auxmat2(2,1)=0.d0
auxmat2(2,2)=hu(n)
auxmat22=matmul(auxmat2,ZZ)
hutilde=1.d0
du1tilde(n-1)=Lambda(1,1)
gu1tilde(n-1,1)=Lambda(1,2)
gu1tilde(n-1,2)=0.d0
du1tilde(N)=Lambda(2,2)
hu1tilde(1,n)=auxmat22(1,1)
hu1tilde(2,n)=auxmat22(2,1)
bu1tilde(1,1,n)=auxmat22(1,2)
bu1tilde(2,1,n)=auxmat22(2,2)
bu1tilde(1,2,n)=0.d0
bu1tilde(2,2,n)=0.d0
hv1tilde(1,n+2)=1.d0
hv1tilde(2,n+2)=0.d0
hv1tilde(3,n+2)=0.d0
hu1tilde(1,n+1)=1.d0
hu1tilde(2,n+1)=0.d0

! *** (8) ****

w1(n-1)=conjg(ZZ(1,1))*gggamma+conjg(ZZ(2,1))*w(n)
w1(n)=conjg(ZZ(1,2))*gggamma+conjg(ZZ(2,2))*w(n)
p1(n-1)=QQ(1,1)*c+QQ(1,2)*p(n)
p1(n)=QQ(2,1)*c+QQ(2,2)*p(n)
q1(n-1)=conjg(ZZ(1,1))*theta+conjg(ZZ(2,1))*q(n)
q1(n)=conjg(ZZ(1,2))*theta+conjg(ZZ(2,2))*q(n)
z1(n)=chi
sigmaanew(n-1)=dv1tilde(n+1)-z1(n)*conjg(w1(n-1))


! *** (9) ****

do j=1,n
   db1(j)=du1tilde(j)-p1(j)*conjg(q1(j))
enddo

! *** (10) ****

do j=1,n
   gv1(j,1)=gv1tilde(j+1,1)
   gv1(j,2)=gv1tilde(j+1,2)
   gv1(j,3)=gv1tilde(j+1,3)
   hv1(1,j)=hv1tilde(1,j+2)
   hv1(2,j)=hv1tilde(2,j+2)
   hv1(3,j)=hv1tilde(3,j+2)
enddo
do j=1,n-1
   bv1(1,1,j)=bv1tilde(1,1,j+2)
   bv1(1,2,j)=bv1tilde(1,2,j+2)
   bv1(1,3,j)=bv1tilde(1,3,j+2)
   bv1(2,1,j)=bv1tilde(2,1,j+2)
   bv1(2,2,j)=bv1tilde(2,2,j+2)
   bv1(2,3,j)=bv1tilde(2,3,j+2)
   bv1(3,1,j)=bv1tilde(3,1,j+2)
   bv1(3,2,j)=bv1tilde(3,2,j+2)
   bv1(3,3,j)=bv1tilde(3,3,j+2)
enddo
do j=2,n
   dv1(j)=dv1tilde(j+1)
enddo
do j=1,n
   gu1(j,1)=gu1tilde(j,1)
   gu1(j,2)=gu1tilde(j,2)
   hu1(1,j)=hu1tilde(1,j+1)
   hu1(2,j)=hu1tilde(2,j+1)
   du1(j)=du1tilde(j)
enddo
do j=2,n-1
   bu1(1,1,j)=bu1tilde(1,1,j+1)
   bu1(1,2,j)=bu1tilde(1,2,j+1)
   bu1(2,1,j)=bu1tilde(2,1,j+1)
   bu1(2,2,j)=bu1tilde(2,2,j+1)
enddo
! Now we have redundant generators for the updated A and B


! ***** CHECK *****
!call recover_A_red(othera,sigmaanew,gv1,hv1,bv1,z1,w1,n)
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(A(i,j)-othera(i,j))*conjg(A(i,j)-othera(i,j))
!      !if (abs(A(i,j)-othera(i,j))>1000*tol) then
!      !   write(*,*)"i,j=",i,j
!      !endif   
!   enddo
!enddo  
!write(*,*)"norm of difference for A_red=",normaF2 
!
!call recover_V_red(othera,sigmaanew,gv1,hv1,bv1,z1,w1,n)
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(V(i,j)-othera(i,j))*conjg(V(i,j)-othera(i,j))
!   enddo
!enddo  
!write(*,*)"norm of difference for V_red=",normaF2 
!
!
! *****************




! *** (11) ****
! Compression of V

Sold(1,1)=gv1(1,1)
Sold(1,2)=gv1(1,2)
Sold(1,3)=gv1(1,3)
gs(1,1)=1.d0
gs(1,2)=0.d0
hs(1,1)=Sold(1,1)*hv1(1,1)+Sold(1,2)*hv1(2,1)+Sold(1,3)*hv1(3,1)
hs(2,1)=0.d0

gs(2,1)=0.d0
gs(2,2)=1.d0
bs(1,1,1)=1.d0
bs(1,2,1)=0.d0
bs(2,1,1)=0.d0
bs(2,2,1)=0.d0
Snew(1,1)=Sold(1,1)*bv1(1,1,1)+Sold(1,2)*bv1(2,1,1)+Sold(1,3)*bv1(3,1,1)
Snew(1,2)=Sold(1,1)*bv1(1,2,1)+Sold(1,2)*bv1(2,2,1)+Sold(1,3)*bv1(3,2,1)
Snew(1,3)=Sold(1,1)*bv1(1,3,1)+Sold(1,2)*bv1(2,3,1)+Sold(1,3)*bv1(3,3,1)
Snew(2,1)=gv1(2,1)
Snew(2,2)=gv1(2,2)
Snew(2,3)=gv1(2,3)
hs(1,2)=Snew(1,1)*hv1(1,2)+Snew(1,2)*hv1(2,2)+Snew(1,3)*hv1(3,2)
hs(2,2)=Snew(2,1)*hv1(1,2)+Snew(2,2)*hv1(2,2)+Snew(2,3)*hv1(3,2)

zvec(1)=gs(1,1)*hs(1,1)+gs(1,2)*hs(2,1)
zvec(2)=dv1(2)
yold=w1(1)


do k=2,n-1
   auxvec(1)=conjg(yold)
   auxvec(2)=conjg(w1(k))
   call givensr(auxvec,VV)
   ynew=auxvec(1)*VV(1,2)+auxvec(2)*VV(2,2)
   auxmat32(1,1)=zvec(1)
   auxmat32(2,1)=zvec(2)
   auxmat32(1,2)=Snew(1,1)*hv1(1,k)+Snew(1,2)*hv1(2,k)+Snew(1,3)*hv1(3,k)
   auxmat32(2,2)=Snew(2,1)*hv1(1,k)+Snew(2,2)*hv1(2,k)+Snew(2,3)*hv1(3,k)
   auxmat32(3,1)=z1(k+1)*conjg(yold)
   auxmat32(3,2)=dv1(k+1)
   AA=matmul(auxmat32,VV)
   yold=ynew
   auxvec3=AA(:,1)
   call hh3(auxvec3,FF)
   Sold=Snew
   
   auxmat(1,1)=Sold(1,1)*bv1(1,1,k)+Sold(1,2)*bv1(2,1,k)+Sold(1,3)*bv1(3,1,k)
   auxmat(1,2)=Sold(1,1)*bv1(1,2,k)+Sold(1,2)*bv1(2,2,k)+Sold(1,3)*bv1(3,2,k)
   auxmat(1,3)=Sold(1,1)*bv1(1,3,k)+Sold(1,2)*bv1(2,3,k)+Sold(1,3)*bv1(3,3,k)
   auxmat(2,1)=Sold(2,1)*bv1(1,1,k)+Sold(2,2)*bv1(2,1,k)+Sold(2,3)*bv1(3,1,k)
   auxmat(2,2)=Sold(2,1)*bv1(1,2,k)+Sold(2,2)*bv1(2,2,k)+Sold(2,3)*bv1(3,2,k)
   auxmat(2,3)=Sold(2,1)*bv1(1,3,k)+Sold(2,2)*bv1(2,3,k)+Sold(2,3)*bv1(3,3,k)
 
   zvec(1)=FF(2,1)*AA(1,2)+FF(2,2)*AA(2,2)+FF(2,3)*AA(3,2)       ! Bb(2,1)
   zvec(2)=FF(3,1)*AA(1,2)+FF(3,2)*AA(2,2)+FF(3,3)*AA(3,2)       ! Bb(3,1)  
   bs(1,1,k)=conjg(FF(2,1))
   bs(1,2,k)=conjg(FF(3,1))
   bs(2,1,k)=conjg(FF(2,2))
   bs(2,2,k)=conjg(FF(3,2))
   gs(k+1,1)=conjg(FF(2,3))
   gs(k+1,2)=conjg(FF(3,3))
   Snew(1,1)=conjg(gs(k+1,1))*gv1(k+1,1)+conjg(bs(1,1,k))*auxmat(1,1)+conjg(bs(2,1,k))*auxmat(2,1)
   Snew(1,2)=conjg(gs(k+1,1))*gv1(k+1,2)+conjg(bs(1,1,k))*auxmat(1,2)+conjg(bs(2,1,k))*auxmat(2,2)
   Snew(1,3)=conjg(gs(k+1,1))*gv1(k+1,3)+conjg(bs(1,1,k))*auxmat(1,3)+conjg(bs(2,1,k))*auxmat(2,3)
   Snew(2,1)=conjg(gs(k+1,2))*gv1(k+1,1)+conjg(bs(1,2,k))*auxmat(1,1)+conjg(bs(2,2,k))*auxmat(2,1)
   Snew(2,2)=conjg(gs(k+1,2))*gv1(k+1,2)+conjg(bs(1,2,k))*auxmat(1,2)+conjg(bs(2,2,k))*auxmat(2,2)
   Snew(2,3)=conjg(gs(k+1,2))*gv1(k+1,3)+conjg(bs(1,2,k))*auxmat(1,3)+conjg(bs(2,2,k))*auxmat(2,3)
   hs(1,k+1)=Snew(1,1)*hv1(1,k+1)+Snew(1,2)*hv1(2,k+1)+Snew(1,3)*hv1(3,k+1)
   hs(2,k+1)=Snew(2,1)*hv1(1,k+1)+Snew(2,2)*hv1(2,k+1)+Snew(2,3)*hv1(3,k+1)

enddo


gvnew=gs
hvnew=hs
bvnew=bs


! *** (12) ****
! Compression of U

xnum=q1(1)
yvec=gu1(1,:)
znum=du1(1)
hunew(1)=0.d0
hunew(2)=gu1(1,1)*hu1(1,2)+gu1(1,2)*hu1(2,2)
gunew(1)=1.d0
do k=2,n-1
  xx(1)=conjg(xnum)
  xx(2)=conjg(q1(k))
  call givensr(xx,WW)
  auxmat2(1,1)=znum
  auxmat2(1,2)=hunew(k)
  auxmat2(2,1)=p1(k)*conjg(xnum)
  auxmat2(2,2)=du1(k)
  Zmat=matmul(auxmat2,WW)
  xnum=xx(1)*WW(1,2)+xx(2)*WW(2,2)
  xx=Zmat(:,1)
  call givens(xx,auxmat22)
 
  bunew(k)=conjg(auxmat22(2,1))
  gunew(k)=conjg(auxmat22(2,2))
  auxvec=yvec
  yvec(1)=conjg(gunew(k))*gu1(k,1)+conjg(bunew(k))*(auxvec(1)*bu1(1,1,k)+auxvec(2)*bu1(2,1,k))
  yvec(2)=conjg(gunew(k))*gu1(k,2)+conjg(bunew(k))*(auxvec(1)*bu1(1,2,k)+auxvec(2)*bu1(2,2,k))
  znum=conjg(gunew(k))*Zmat(2,2)+conjg(bunew(k))*Zmat(1,2)
  hunew(k+1)=yvec(1)*hu1(1,k+1)+yvec(2)*hu1(2,k+1)
enddo

! update
sigmaa=sigmaanew
gv=gvnew
hv=hvnew
bv=bvnew
gu=gunew
hu=hunew
bu=bunew
db=db1
p=p1
q=q1
z=z1
w=w1

! ***** CHECK *****
!call recover_A(othera,sigmaanew,gvnew,hvnew,bvnew,z1,w1,n)
!call recover_B(otherb,db1,gunew,hunew,bunew,p1,q1,n)
!write(*,*)"A=",A
!write(*,*)"othera=",othera
!write(*,*)"B=",B
!write(*,*)"otherb=",otherb 
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(A(i,j)-othera(i,j))*conjg(A(i,j)-othera(i,j))
!   enddo
!enddo  
!normbF2=0.d0
!do i=1,n
!   do j=1,n
!      normbF2=normbF2+(B(i,j)-otherb(i,j))*conjg(B(i,j)-otherb(i,j))
!   enddo
!enddo  
!write(*,*)"norm of difference for A=",normaF2 
!write(*,*)"norm of difference for B=",normbF2 
! *****************



!**********************************
! if kfin is small, stop loop early
! (update of generators will be incomplete,
! checks on whole matrices won't be useful here)
else
!**********************************

do j=1,n-1
   sigmav(j)=sigmaa(j)+z(j+1)*conjg(w(j))
enddo   
epsi=gv(1,1)*hv(1,1)+gv(1,2)*hv(2,1)-z(1)*conjg(w(1))
if (kin==1) then
    QQ=Qin
else
    QQ(1,1)=1.d0
    QQ(2,2)=1.d0
    QQ(1,2)=0.d0
    QQ(2,1)=0.d0
endif

! *** (1) ****

GGamma=0.d0
auxmat(1,1)=gv(1,1)*hv(1,1)+gv(1,2)*hv(2,1)
auxmat(1,2)=gv(1,1)*bv(1,1,1)+gv(1,2)*bv(2,1,1)
auxmat(1,3)=gv(1,1)*bv(1,2,1)+gv(1,2)*bv(2,2,1)
auxmat(2,1)=sigmav(1)
auxmat(2,2)=gv(2,1)
auxmat(2,3)=gv(2,2)
GGamma(:,1:3)=matmul(QQ,auxmat)

gv1tilde(2,:)=GGamma(1,1:3)
fv=GGamma(2,1)
phiv=GGamma(2,2:3)
auxvec(1)=QQ(1,1)*z(1)+QQ(1,2)*z(2)
auxvec(2)=QQ(2,1)*z(1)+QQ(2,2)*z(2)

z1(1)=auxvec(1)
chi=auxvec(2)
fa=fv-chi*conjg(w(1))
phia=phiv

! ******* CHECK ******************
!mat2n=A(1:2,:)
!A(1:2,:)=matmul(QQ,mat2n)
!mat2n=V(1:2,:)
!V(1:2,:)=matmul(QQ,mat2n)
! ********************************


! *** (2) ****

gggamma=w(1)
c=p(1)
theta=q(1)
du(1)=db(1)+p(1)*conjg(q(1))
fu=du(1)
phiu=gu(1)
fb=db(1)
phi=gu(1)

! *** (3) ****

do k=1,kfin+3
   ! *** (3a) ****
   vareps=phi*hu(k+1)-c*conjg(q(k+1))
   epsi=phia(1)*hv(1,k+1)+phia(2)*hv(2,k+1)-chi*conjg(w(k+1))
   du(k+1)=db(k+1)+p(k+1)*conjg(q(k+1))
   auxmat2(1,1)=fb
   auxmat2(1,2)=vareps
   auxmat2(2,1)=0.d0
   auxmat2(2,2)=db(k+1)
   PPhi=matmul(QQ,auxmat2)
   
   
   ! ******* CHECK ******************
   !mat2n=B(k:k+1,:)
   !B(k:k+1,:)=matmul(QQ,mat2n)
   !mat2n=U(k:k+1,:)
   !U(k:k+1,:)=matmul(QQ,mat2n)
   ! ********************************
   
   
   ! *** (3b) ****
   if (k>=kin) then
      xx(1)=conjg(PPhi(2,1))
      xx(2)=conjg(PPhi(2,2))
      call givens(xx,GG)
      ZZ(1,1)=conjg(GG(2,1))
      ZZ(1,2)=conjg(GG(1,1))
      ZZ(2,1)=conjg(GG(2,2))
      ZZ(2,2)=conjg(GG(1,2))
   else
      ZZ(1,1)=1.d0
      ZZ(1,2)=0.d0
      ZZ(2,1)=0.d0
      ZZ(2,2)=1.d0
   endif
   auxmat2(1,1)=fa
   auxmat2(1,2)=epsi
   auxmat2(2,1)=0.d0
   auxmat2(2,2)=sigmaa(k+1)
   
   Omega=matmul(auxmat2,ZZ)   
   
   ! ****** CHECK ********
   !matn2=A(:,k:k+1)
   !A(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=V(:,k:k+1)
   !V(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=B(:,k:k+1)
   !B(:,k:k+1)=matmul(matn2,ZZ)
   !matn2=U(:,k:k+1)
   !U(:,k:k+1)=matmul(matn2,ZZ)
   ! **********************
   
   
   ! *** (3c) ****
   if (k>=kin .AND. k<=kfin-2) then
      xx=Omega(:,1)
      call givens(xx,newQ)
      sigmaanew(k)=newQ(1,1)*xx(1)+newQ(1,2)*xx(2)
   else
      if (k+1==kin) then
         newQ=Qin
      else
         newQ(1,1)=1.d0
         newQ(1,2)=0.d0
         newQ(2,1)=0.d0
         newQ(2,2)=1.d0                     
      end if
      sigmaanew(k)=newQ(1,1)*Omega(1,1)+newQ(1,2)*Omega(2,1)
   endif
   
   
   ! ***** CHECK *********
   !mat2n=A(k+1:k+2,:)
   !A(k+1:k+2,:)=matmul(newQ,mat2n)
   !mat2n=V(k+1:k+2,:)
   !V(k+1:k+2,:)=matmul(newQ,mat2n)
   ! *********************
   
   
   
   ! *** (3d) ****
   auxmat24(1,1)=fv
   auxmat24(1,2)=phiv(1)*hv(1,k+1)+phiv(2)*hv(2,k+1)
   auxmat24(1,3)=phiv(1)*bv(1,1,k+1)+phiv(2)*bv(2,1,k+1)
   auxmat24(1,4)=phiv(1)*bv(1,2,k+1)+phiv(2)*bv(2,2,k+1)
   auxmat24(2,1)=z(k+2)*conjg(gggamma)
   auxmat24(2,2)=sigmav(k+1)
   auxmat24(2,3)=gv(k+2,1)
   auxmat24(2,4)=gv(k+2,2)
   Gammap=matmul(newQ,auxmat24)
   GGamma=Gammap
   
   auxmat2(1,1)=Gammap(1,1)*ZZ(1,1)+Gammap(1,2)*ZZ(2,1)
   auxmat2(1,2)=Gammap(1,1)*ZZ(1,2)+Gammap(1,2)*ZZ(2,2)
   auxmat2(2,1)=Gammap(2,1)*ZZ(1,1)+Gammap(2,2)*ZZ(2,1)
   auxmat2(2,2)=Gammap(2,1)*ZZ(1,2)+Gammap(2,2)*ZZ(2,2)
   
   GGamma(1,1)=auxmat2(1,1)
   GGamma(1,2)=auxmat2(1,2)
   GGamma(2,1)=auxmat2(2,1)
   GGamma(2,2)=auxmat2(2,2)
   auxmat34(1,1)=1.d0
   auxmat34(1,2)=0.d0
   auxmat34(1,3)=0.d0
   auxmat34(1,4)=0.d0
   auxmat34(2,1)=0.d0
   auxmat34(3,1)=0.d0
   auxmat34(2:3,2)=hv(:,k+1)
   auxmat34(2:3,3:4)=bv(:,:,k+1)
   CC=auxmat34
   auxmat32=matmul(auxmat34(:,1:2),ZZ)
   CC(1,1)=auxmat32(1,1)
   CC(1,2)=auxmat32(1,2)
   CC(2,1)=auxmat32(2,1)
   CC(2,2)=auxmat32(2,2)
   CC(3,1)=auxmat32(3,1)
   CC(3,2)=auxmat32(3,2)
   dv1tilde(k+2)=GGamma(1,1)
   gv1tilde(k+2,:)=GGamma(1,2:4)
   betav=GGamma(2,2:4)
   fv=betav(1)
   phiv=betav(2:3)
   hv1tilde(1,k+2)=CC(1,1)
   hv1tilde(2,k+2)=CC(2,1)
   hv1tilde(3,k+2)=CC(3,1)
   bv1tilde(1,1,k+2)=CC(1,2)
   bv1tilde(1,2,k+2)=CC(1,3)
   bv1tilde(1,3,k+2)=CC(1,4)
   bv1tilde(2,1,k+2)=CC(2,2)
   bv1tilde(2,2,k+2)=CC(2,3)
   bv1tilde(2,3,k+2)=CC(2,4)
   bv1tilde(3,1,k+2)=CC(3,2)
   bv1tilde(3,2,k+2)=CC(3,3)
   bv1tilde(3,3,k+2)=CC(3,4)
   ! *** (3e) ****
   auxmat(1,1)=fu
   auxmat(1,2)=phiu*hu(k+1)
   auxmat(1,3)=phiu*bu(k+1)
   auxmat(2,1)=p(k+1)*conjg(theta)
   auxmat(2,2)=du(k+1)
   auxmat(2,3)=gu(k+1)
   auxmat2=matmul(auxmat(:,1:2),ZZ)
   auxmat(1,1)=auxmat2(1,1)
   auxmat(1,2)=auxmat2(1,2)
   auxmat(2,1)=auxmat2(2,1)
   auxmat(2,2)=auxmat2(2,2)
   Lambda=matmul(QQ,auxmat)
   auxmat(1,1)=1.d0
   auxmat(1,2)=0.d0
   auxmat(1,3)=0.d0
   auxmat(2,1)=0.d0
   auxmat(2,2)=hu(k+1)
   auxmat(2,3)=bu(k+1)
   auxmat2=matmul(auxmat(:,1:2),ZZ)
   auxmat(1,1)=auxmat2(1,1)
   auxmat(1,2)=auxmat2(1,2)
   auxmat(2,1)=auxmat2(2,1)
   auxmat(2,2)=auxmat2(2,2)
   du1tilde(k)=Lambda(1,1)
   gu1tilde(k,1)=Lambda(1,2)
   gu1tilde(k,2)=Lambda(1,3)
   betau=Lambda(2,2:3)
   fu=betau(1)
   phiu=betau(2)
   hu1tilde(1,k+1)=auxmat(1,1)
   hu1tilde(2,k+1)=auxmat(2,1)
   bu1tilde(1,1,k+1)=auxmat(1,2)
   bu1tilde(1,2,k+1)=auxmat(1,3)
   bu1tilde(2,1,k+1)=auxmat(2,2)
   bu1tilde(2,2,k+1)=auxmat(2,3)
   ! *** (3f) ****
   auxvec(1)=newQ(1,1)*chi+newQ(1,2)*z(k+2)
   auxvec(2)=newQ(2,1)*chi+newQ(2,2)*z(k+2)
   z1(k+1)=auxvec(1)
   chi=auxvec(2)
   auxvec(1)=conjg(ZZ(1,1))*gggamma+conjg(ZZ(2,1))*w(k+1)
   auxvec(2)=conjg(ZZ(1,2))*gggamma+conjg(ZZ(2,2))*w(k+1)   
   w1(k)=auxvec(1)
   gggamma=auxvec(2)
   ! *** (3g) ****
   auxvec(1)=conjg(ZZ(1,1))*theta+conjg(ZZ(2,1))*q(k+1)
   auxvec(2)=conjg(ZZ(1,2))*theta+conjg(ZZ(2,2))*q(k+1)
   q1(k)=auxvec(1)
   theta=auxvec(2)
   auxvec(1)=QQ(1,1)*c+QQ(1,2)*p(k+1)
   auxvec(2)=QQ(2,1)*c+QQ(2,2)*p(k+1)
   p1(k)=auxvec(1)
   c=auxvec(2)
   ! *** (3h) ****
   fa=fv-chi*conjg(gggamma)
   phia=phiv
   fb=fu-c*conjg(theta)
   phi=phiu
   QQ=newQ  
   
enddo

   

! *** (4) ****

vareps=phi*hu(n)-c*conjg(q(n))
du(n)=db(n)+p(n)*conjg(q(n))
auxmat2(1,1)=fb
auxmat2(1,2)=vareps
auxmat2(2,1)=0.d0
auxmat2(2,2)=db(n)
PPhi=matmul(QQ,auxmat2)

! ***** CHECK ******
!mat2n=B(n-1:n,:)
!B(n-1:n,:)=matmul(QQ,mat2n)
!mat2n=U(n-1:n,:)
!U(n-1:n,:)=matmul(QQ,mat2n)
! ******************


! *** (5) ****

xx(1)=conjg(PPhi(2,1))
xx(2)=conjg(PPhi(2,2))
call givens(xx,GG)
ZZ(1,1)=conjg(GG(2,1))
ZZ(1,2)=conjg(GG(1,1))
ZZ(2,1)=conjg(GG(2,2))
ZZ(2,2)=conjg(GG(1,2))

! *** (6) ****
! empty, no need for last updates


! ***** CHECK *****
!matn2=A(:,n-1:n)
!A(:,n-1:n)=matmul(matn2,ZZ)
!matn2=V(:,n-1:n)
!V(:,n-1:n)=matmul(matn2,ZZ)
!matn2=B(:,n-1:n)
!B(:,n-1:n)=matmul(matn2,ZZ)
!matn2=U(:,n-1:n)
!U(:,n-1:n)=matmul(matn2,ZZ)
! *****************




! *** (7) ****
! empty, no need for last updates

! *** (8) ****
! empty, no need for last updates

! *** (9) ****

do j=1,kfin+3
   db1(j)=du1tilde(j)-p1(j)*conjg(q1(j))
enddo

! *** (10) ****

do j=1,kfin+3
   gv1(j,1)=gv1tilde(j+1,1)
   gv1(j,2)=gv1tilde(j+1,2)
   gv1(j,3)=gv1tilde(j+1,3)
   hv1(1,j)=hv1tilde(1,j+2)
   hv1(2,j)=hv1tilde(2,j+2)
   hv1(3,j)=hv1tilde(3,j+2)
enddo
do j=1,kfin+3
   bv1(1,1,j)=bv1tilde(1,1,j+2)
   bv1(1,2,j)=bv1tilde(1,2,j+2)
   bv1(1,3,j)=bv1tilde(1,3,j+2)
   bv1(2,1,j)=bv1tilde(2,1,j+2)
   bv1(2,2,j)=bv1tilde(2,2,j+2)
   bv1(2,3,j)=bv1tilde(2,3,j+2)
   bv1(3,1,j)=bv1tilde(3,1,j+2)
   bv1(3,2,j)=bv1tilde(3,2,j+2)
   bv1(3,3,j)=bv1tilde(3,3,j+2)
enddo
do j=2,kfin+3
   dv1(j)=dv1tilde(j+1)
enddo
do j=1,kfin+3
   gu1(j,1)=gu1tilde(j,1)
   gu1(j,2)=gu1tilde(j,2)
   hu1(1,j)=hu1tilde(1,j+1)
   hu1(2,j)=hu1tilde(2,j+1)
   du1(j)=du1tilde(j)
enddo
do j=2,kfin+3
   bu1(1,1,j)=bu1tilde(1,1,j+1)
   bu1(1,2,j)=bu1tilde(1,2,j+1)
   bu1(2,1,j)=bu1tilde(2,1,j+1)
   bu1(2,2,j)=bu1tilde(2,2,j+1)
enddo
! Now we have redundant generators for the updated A and B


! ***** CHECK *****
!call recover_A_red(othera,sigmaanew,gv1,hv1,bv1,z1,w1,n)
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(A(i,j)-othera(i,j))*conjg(A(i,j)-othera(i,j))
!      !if (abs(A(i,j)-othera(i,j))>1000*tol) then
!      !   write(*,*)"i,j=",i,j
!      !endif   
!   enddo
!enddo  
!write(*,*)"norm of difference for A_red=",normaF2 
!
!call recover_V_red(othera,sigmaanew,gv1,hv1,bv1,z1,w1,n)
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(V(i,j)-othera(i,j))*conjg(V(i,j)-othera(i,j))
!   enddo
!enddo  
!write(*,*)"norm of difference for V_red=",normaF2 
!
!
! *****************




! *** (11) ****
! Compression of V

Sold(1,1)=gv1(1,1)
Sold(1,2)=gv1(1,2)
Sold(1,3)=gv1(1,3)
gs(1,1)=1.d0
gs(1,2)=0.d0
hs(1,1)=Sold(1,1)*hv1(1,1)+Sold(1,2)*hv1(2,1)+Sold(1,3)*hv1(3,1)
hs(2,1)=0.d0

gs(2,1)=0.d0
gs(2,2)=1.d0
bs(1,1,1)=1.d0
bs(1,2,1)=0.d0
bs(2,1,1)=0.d0
bs(2,2,1)=0.d0
Snew(1,1)=Sold(1,1)*bv1(1,1,1)+Sold(1,2)*bv1(2,1,1)+Sold(1,3)*bv1(3,1,1)
Snew(1,2)=Sold(1,1)*bv1(1,2,1)+Sold(1,2)*bv1(2,2,1)+Sold(1,3)*bv1(3,2,1)
Snew(1,3)=Sold(1,1)*bv1(1,3,1)+Sold(1,2)*bv1(2,3,1)+Sold(1,3)*bv1(3,3,1)
Snew(2,1)=gv1(2,1)
Snew(2,2)=gv1(2,2)
Snew(2,3)=gv1(2,3)
hs(1,2)=Snew(1,1)*hv1(1,2)+Snew(1,2)*hv1(2,2)+Snew(1,3)*hv1(3,2)
hs(2,2)=Snew(2,1)*hv1(1,2)+Snew(2,2)*hv1(2,2)+Snew(2,3)*hv1(3,2)

zvec(1)=gs(1,1)*hs(1,1)+gs(1,2)*hs(2,1)
zvec(2)=dv1(2)
yold=w1(1)


do k=2,kfin+3
   auxvec(1)=conjg(yold)
   auxvec(2)=conjg(w1(k))
   call givensr(auxvec,VV)
  
   ynew=auxvec(1)*VV(1,2)+auxvec(2)*VV(2,2)
   auxmat32(1,1)=zvec(1)
   auxmat32(2,1)=zvec(2)
   auxmat32(1,2)=Snew(1,1)*hv1(1,k)+Snew(1,2)*hv1(2,k)+Snew(1,3)*hv1(3,k)
   auxmat32(2,2)=Snew(2,1)*hv1(1,k)+Snew(2,2)*hv1(2,k)+Snew(2,3)*hv1(3,k)
   auxmat32(3,1)=z1(k+1)*conjg(yold)
   auxmat32(3,2)=dv1(k+1)
   AA=matmul(auxmat32,VV)
   yold=ynew
   auxvec3=AA(:,1)
   call hh3(auxvec3,FF)
   Sold=Snew
   
   auxmat(1,1)=Sold(1,1)*bv1(1,1,k)+Sold(1,2)*bv1(2,1,k)+Sold(1,3)*bv1(3,1,k)
   auxmat(1,2)=Sold(1,1)*bv1(1,2,k)+Sold(1,2)*bv1(2,2,k)+Sold(1,3)*bv1(3,2,k)
   auxmat(1,3)=Sold(1,1)*bv1(1,3,k)+Sold(1,2)*bv1(2,3,k)+Sold(1,3)*bv1(3,3,k)
   auxmat(2,1)=Sold(2,1)*bv1(1,1,k)+Sold(2,2)*bv1(2,1,k)+Sold(2,3)*bv1(3,1,k)
   auxmat(2,2)=Sold(2,1)*bv1(1,2,k)+Sold(2,2)*bv1(2,2,k)+Sold(2,3)*bv1(3,2,k)
   auxmat(2,3)=Sold(2,1)*bv1(1,3,k)+Sold(2,2)*bv1(2,3,k)+Sold(2,3)*bv1(3,3,k)
 
   zvec(1)=FF(2,1)*AA(1,2)+FF(2,2)*AA(2,2)+FF(2,3)*AA(3,2)       ! Bb(2,1)
   zvec(2)=FF(3,1)*AA(1,2)+FF(3,2)*AA(2,2)+FF(3,3)*AA(3,2)       ! Bb(3,1)  
   bs(1,1,k)=conjg(FF(2,1))
   bs(1,2,k)=conjg(FF(3,1))
   bs(2,1,k)=conjg(FF(2,2))
   bs(2,2,k)=conjg(FF(3,2))
   gs(k+1,1)=conjg(FF(2,3))
   gs(k+1,2)=conjg(FF(3,3))
   Snew(1,1)=conjg(gs(k+1,1))*gv1(k+1,1)+conjg(bs(1,1,k))*auxmat(1,1)+conjg(bs(2,1,k))*auxmat(2,1)
   Snew(1,2)=conjg(gs(k+1,1))*gv1(k+1,2)+conjg(bs(1,1,k))*auxmat(1,2)+conjg(bs(2,1,k))*auxmat(2,2)
   Snew(1,3)=conjg(gs(k+1,1))*gv1(k+1,3)+conjg(bs(1,1,k))*auxmat(1,3)+conjg(bs(2,1,k))*auxmat(2,3)
   Snew(2,1)=conjg(gs(k+1,2))*gv1(k+1,1)+conjg(bs(1,2,k))*auxmat(1,1)+conjg(bs(2,2,k))*auxmat(2,1)
   Snew(2,2)=conjg(gs(k+1,2))*gv1(k+1,2)+conjg(bs(1,2,k))*auxmat(1,2)+conjg(bs(2,2,k))*auxmat(2,2)
   Snew(2,3)=conjg(gs(k+1,2))*gv1(k+1,3)+conjg(bs(1,2,k))*auxmat(1,3)+conjg(bs(2,2,k))*auxmat(2,3)
   hs(1,k+1)=Snew(1,1)*hv1(1,k+1)+Snew(1,2)*hv1(2,k+1)+Snew(1,3)*hv1(3,k+1)
   hs(2,k+1)=Snew(2,1)*hv1(1,k+1)+Snew(2,2)*hv1(2,k+1)+Snew(2,3)*hv1(3,k+1)

enddo


gvnew=gs
hvnew=hs
bvnew=bs


! *** (12) ****
! Compression of U

xnum=q1(1)
yvec=gu1(1,:)
znum=du1(1)
hunew(1)=0.d0
hunew(2)=gu1(1,1)*hu1(1,2)+gu1(1,2)*hu1(2,2)
gunew(1)=1.d0
do k=2,kfin+3
  xx(1)=conjg(xnum)
  xx(2)=conjg(q1(k))
  call givensr(xx,WW)
  auxmat2(1,1)=znum
  auxmat2(1,2)=hunew(k)
  auxmat2(2,1)=p1(k)*conjg(xnum)
  auxmat2(2,2)=du1(k)
  Zmat=matmul(auxmat2,WW)
  xnum=xx(1)*WW(1,2)+xx(2)*WW(2,2)
  xx=Zmat(:,1)
  call givens(xx,auxmat22)
  
  bunew(k)=conjg(auxmat22(2,1))
  gunew(k)=conjg(auxmat22(2,2))
  auxvec=yvec
  yvec(1)=conjg(gunew(k))*gu1(k,1)+conjg(bunew(k))*(auxvec(1)*bu1(1,1,k)+auxvec(2)*bu1(2,1,k))
  yvec(2)=conjg(gunew(k))*gu1(k,2)+conjg(bunew(k))*(auxvec(1)*bu1(1,2,k)+auxvec(2)*bu1(2,2,k))
  znum=conjg(gunew(k))*Zmat(2,2)+conjg(bunew(k))*Zmat(1,2)
  hunew(k+1)=yvec(1)*hu1(1,k+1)+yvec(2)*hu1(2,k+1)
enddo

! update
sigmaa=sigmaanew
gv=gvnew
hv=hvnew
bv=bvnew
gu=gunew
hu=hunew
bu=bunew
db=db1
p=p1
q=q1
z=z1
w=w1

! ***** CHECK *****
!call recover_A(othera,sigmaanew,gvnew,hvnew,bvnew,z1,w1,n)
!call recover_B(otherb,db1,gunew,hunew,bunew,p1,q1,n)
!write(*,*)"A=",A
!write(*,*)"othera=",othera
!write(*,*)"B=",B
!write(*,*)"otherb=",otherb 
!normaF2=0.d0
!do i=1,n
!   do j=1,n
!      normaF2=normaF2+(A(i,j)-othera(i,j))*conjg(A(i,j)-othera(i,j))
!   enddo
!enddo  
!normbF2=0.d0
!do i=1,n
!   do j=1,n
!      normbF2=normbF2+(B(i,j)-otherb(i,j))*conjg(B(i,j)-otherb(i,j))
!   enddo
!enddo  
!write(*,*)"norm of difference for A=",normaF2 
!write(*,*)"norm of difference for B=",normbF2 
! *****************


!***************
endif
!***************


!*************************************
! end of structured QZ sweep
!*************************************




      numit=numit+1
      total_it=total_it+1
      !write(*,*) "."
   end if
enddo
adiag(1)=gv(1,1)*hv(1,1)+gv(1,2)*hv(2,1)-z(1)*conjg(w(1))
bdiag(1)=db(1)

do j=1,n
   myeig(j)=adiag(j)/bdiag(j)
enddo

return

end subroutine fast_companion_QZ


!************************************************************


subroutine recover_A(A,sigmaa,gv,hv,bv,z,w,n)

! Recovers the matrix A from its quasiseparable generators.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: sigmaa(n-1),gv(n,2),hv(2,n),bv(2,2,n-1),z(n),w(n)
complex(kind(1d0)), intent(out) :: A(n,n)

integer :: j,k
complex(kind(1d0)) :: x(2),auxx(2)

A=0.d0

! Fill subdiagonal:
do k=1,n-1
   A(k+1,k)=sigmaa(k)
enddo

! Fill upper triangular part:
do j=1,n-1
   x=gv(j,:)
   do k=j,n-1
      A(j,k)=x(1)*hv(1,k)+x(2)*hv(2,k)-z(j)*conjg(w(k))
      auxx=x
      x(1)=auxx(1)*bv(1,1,k)+auxx(2)*bv(2,1,k)
      x(2)=auxx(1)*bv(1,2,k)+auxx(2)*bv(2,2,k)
   enddo
   A(j,n)=x(1)*hv(1,n)+x(2)*hv(2,n)-z(j)*conjg(w(n))
enddo
A(n,n)=gv(n,1)*hv(1,n)+gv(n,2)*hv(2,n)-z(n)*conjg(w(n))

return

end subroutine recover_A


!************************************************************


subroutine recover_A_red(A,sigmaa,gv3,hv3,bv3,z,w,n)

! Recovers the matrix A from its redundant quasiseparable generators.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: sigmaa(n-1),gv3(n,3),hv3(3,n),bv3(3,3,n),z(n),w(n)
complex(kind(1d0)), intent(out) :: A(n,n)

integer :: j,k
complex(kind(1d0)) :: x(3),auxx(3)

A=0.d0

! Fill subdiagonal:
do k=1,n-1
   A(k+1,k)=sigmaa(k)
enddo

! Fill upper triangular part:
do j=1,n-1
   x=gv3(j,:)
   do k=j,n-1
      A(j,k)=x(1)*hv3(1,k)+x(2)*hv3(2,k)+x(3)*hv3(3,k)-z(j)*conjg(w(k))
      auxx=x
      x(1)=auxx(1)*bv3(1,1,k)+auxx(2)*bv3(2,1,k)+auxx(3)*bv3(3,1,k)
      x(2)=auxx(1)*bv3(1,2,k)+auxx(2)*bv3(2,2,k)+auxx(3)*bv3(3,2,k)
      x(3)=auxx(1)*bv3(1,3,k)+auxx(2)*bv3(2,3,k)+auxx(3)*bv3(3,3,k)
   enddo
   A(j,n)=x(1)*hv3(1,n)+x(2)*hv3(2,n)+x(3)*hv3(3,n)-z(j)*conjg(w(n))
enddo
A(n,n)=gv3(n,1)*hv3(1,n)+gv3(n,2)*hv3(2,n)+gv3(n,3)*hv3(3,n)-z(n)*conjg(w(n))

return

end subroutine recover_A_red


!************************************************************


subroutine recover_V_red(V,sigmaa,gv3,hv3,bv3,z,w,n)

! Recovers the matrix V from its redundant quasiseparable generators.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: sigmaa(n-1),gv3(n,3),hv3(3,n),bv3(3,3,n),z(n),w(n)
complex(kind(1d0)), intent(out) :: V(n,n)

integer :: j,k
complex(kind(1d0)) :: x(3),auxx(3)

V=0.d0

! Fill subdiagonal:
do k=1,n-1
   V(k+1,k)=sigmaa(k)+z(k+1)*conjg(w(k))
enddo

! Fill upper triangular part:
do j=1,n-1
   x=gv3(j,:)
   do k=j,n-1
      V(j,k)=x(1)*hv3(1,k)+x(2)*hv3(2,k)+x(3)*hv3(3,k)
      auxx=x
      x(1)=auxx(1)*bv3(1,1,k)+auxx(2)*bv3(2,1,k)+auxx(3)*bv3(3,1,k)
      x(2)=auxx(1)*bv3(1,2,k)+auxx(2)*bv3(2,2,k)+auxx(3)*bv3(3,2,k)
      x(3)=auxx(1)*bv3(1,3,k)+auxx(2)*bv3(2,3,k)+auxx(3)*bv3(3,3,k)
   enddo
   V(j,n)=x(1)*hv3(1,n)+x(2)*hv3(2,n)+x(3)*hv3(3,n)
enddo
V(n,n)=gv3(n,1)*hv3(1,n)+gv3(n,2)*hv3(2,n)+gv3(n,3)*hv3(3,n)

! Fill lower triangular part:
do j=3,n
   do k=1,j-2
      V(j,k)=z(j)*conjg(w(k))
   enddo
enddo

return

end subroutine recover_V_red


!************************************************************

subroutine recover_B(B,db,gu,hu,bu,p,q,n)

! Recovers the matrix B from its quasiseparable generators.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: db(n),gu(n-1),hu(n),bu(n-1),p(n),q(n)
complex(kind(1d0)), intent(out) :: B(n,n)

integer :: j,k
complex(kind(1d0)) :: x

B=0.d0

! Fill diagonal:
do j=1,n
   B(j,j)=db(j)
enddo

! Fill upper triangular part:
do j=1,n-2
   x=gu(j)
   do k=j+1,n-1
      B(j,k)=x*hu(k)-p(j)*conjg(q(k))
      x=x*bu(k)
   enddo
   B(j,n)=x*hu(n)-p(j)*conjg(q(n))
enddo
B(n-1,n)=gu(n-1)*hu(n)-p(n-1)*conjg(q(n))
  

return

end subroutine

!************************************************************



subroutine hh3(x,H)

! Computes 3x3 Householder matrix that reduces the vector x.

implicit none

complex(kind(1d0)), intent(in) :: x(3)
complex(kind(1d0)), intent(out) :: H(3,3)
complex(kind(1d0)) :: xx(2),yy(2),GG1(2,2),GG2(2,2)

xx(1)=x(2)
xx(2)=x(3)
call givens(xx,GG1)
yy(1)=x(1)
yy(2)=GG1(1,1)*xx(1)+GG1(1,2)*xx(2)
call givens(yy,GG2)
H(1,1)=GG2(1,1)
H(1,2)=GG2(1,2)*GG1(1,1)
H(1,3)=GG2(1,2)*GG1(1,2)
H(2,1)=GG2(2,1)
H(2,2)=GG2(2,2)*GG1(1,1)
H(2,3)=GG2(2,2)*GG1(1,2)
H(3,1)=0.d0
H(3,2)=GG1(2,1)
H(3,3)=GG1(2,2)

return

end subroutine hh3



!************************************************************


subroutine companion_generators(n,poly,sigmaa,gv,hv,bv,z,w,db,gu,hu,bu,p,q)

! Computes quasiseparable generators for the companion pencil
! associated with the polynomial poly.

implicit none

integer, intent(in) :: n
complex(kind(1d0)), intent(in) :: poly(n+1)
complex(kind(1d0)), intent(out) :: sigmaa(n-1),gv(n,2),hv(2,n),bv(2,2,n-1),z(n),w(n)
complex(kind(1d0)), intent(out) :: db(n),gu(n-1),hu(n),bu(n-1),p(n),q(n)
integer :: j

 
! Compute generators for A. 
do j=1,n-1
   sigmaa(j)=1.d0
   bv(1,1,j)=1.d0
   bv(1,2,j)=0.d0
   bv(2,1,j)=0.d0
   bv(2,2,j)=1.d0
enddo
gv=0.d0
gv(1,1)=1.d0
hv=0.d0
hv(1,n)=1.d0
z=0.d0
z(1)=1.d0
w=conjg(poly(2:n+1))
w(n)=w(n)+1

! Compute generators for B.
db=1.d0
db(1)=poly(1)
gu=0.d0
hu=0.d0
bu=0.d0
p=0.d0
p(1)=1.d0
q=0.d0
q(1)=1.d0-conjg(poly(1))

return

end subroutine companion_generators


!********************************************************

subroutine givens(x,G)

! Computes 2x2 Givens matrix such that the second
! entry of G*x is zero.

implicit none
complex(kind(1d0)), intent(in) :: x(2)
complex(kind(1d0)), intent(out) :: G(2,2)
complex(kind(1d0)) :: w1,w2
real(kind(1d0)) :: s

w1=x(1)
w2=x(2)

if (w2==0.d0) then
G(1,1)=1.d0
G(2,2)=1.d0
G(1,2)=0.d0
G(2,1)=0.d0
else

s=sqrt(w1*conjg(w1)+w2*conjg(w2))
G(1,1)=conjg(w1)/s
G(1,2)=conjg(w2)/s
G(2,1)=-w2/s
G(2,2)=w1/s

endif

return

end subroutine givens 


!**********************************************************

subroutine givensr(y,H)

! Computes 2x2 Givens matrix that acts on the right
! of row vector y and puts a zero in the leftmost entry.

implicit none
complex(kind(1d0)), intent(in) :: y(2)
complex(kind(1d0)), intent(out) :: H(2,2)
complex(kind(1d0)) :: G(2,2) 

call givens(y,G)
H(1,1)=G(2,1)
H(1,2)=G(1,1)
H(2,1)=G(2,2)
H(2,2)=G(1,2)

return

end subroutine givensr


!*********************************************************



subroutine	sort_eig(n,eR,eI,seR,seI)
  !!
  !!	sort the sequence according to real part first, if two complex numbers
  !!		have the same real part, put the one with bigger imaginary part in
  !!		front
  !!

  implicit none
  integer				n
  double precision	eR(n),eI(n),seR(n),seI(n),tmpR,tmpI
  integer				i,IY(n),KFLAG,tmpindx

  do i = 1,n
     IY(i) = i
  enddo

  seR = eR
  call ssort(seR,iy,n,kflag)

  do i = 1,n
     seI(i) = eI(iy(i))
  enddo

  tmpR = seR(1)
  tmpI = seI(1)
  do i = 2,n
     if (tmpR == seR(i)) then
        if (tmpI < seI(i)) then
           seI(i-1) = seI(i)
           seI(i)	= tmpI
           tmpindx	= IY(i-1)
           IY(i-1)	= IY(i)
           IY(i) = tmpindx
        endif
     else
        tmpR = seR(i)
        tmpI = seI(i)
     endif
  enddo

  return

end subroutine sort_eig




!**********************************************************


SUBROUTINE init_random_seed()
            INTEGER :: i, n, clock
            INTEGER, DIMENSION(:), ALLOCATABLE :: seed
          
            CALL RANDOM_SEED(size = n)
            ALLOCATE(seed(n))
          
            CALL SYSTEM_CLOCK(COUNT=clock)
          
            seed = clock + 37 * (/ (i - 1, i = 1, n) /)
            CALL RANDOM_SEED(PUT = seed)
          
            DEALLOCATE(seed)
          END SUBROUTINE
