module genschur

! Written by Paul Klein in May 2002
!
! Calls: zgges from LAPACK, see http://www.netlib.org/lapack

contains

subroutine zsolab(f,p,a,b,n,k,retco)

implicit none

integer i,j,k,l,m,n
real*8, dimension(n,n) :: a,b
complex*16, dimension(n,n) :: s,t
complex*16 alpha(n),beta(n),vsl(n,n),vsr(n,n)
complex*16,allocatable :: work(:)
complex*16, dimension(k,k) :: z11, z11i
complex*16, dimension(k,k) :: s11i, t11
complex*16 z21(n-k,k)
real*8 :: f(n-k,k),p(k,k)
integer retco

character*1 jobvsl,jobvsr,sort
integer lda,ldb,sdim,ldvsl,ldvsr,lwork,rwork,info
logical bwork(n)

! Executable section starts here

retco = 0  ! Everything is fine so far
lwork = -1
rwork = 8*n
allocate(work(1))
s = dcmplx(a)
t = dcmplx(b)

ldvsl = n
ldvsr = n
jobvsl = 'n'  ! No, do not calculate left Schur vectors (Q)
jobvsr = 'v'  ! Yes, calculate right Schur vectors (Z)
sort = 's'    ! Yes, order the eigenvalues
lda = n
ldb = n
ldvsl = n
ldvsr = n


! Find the optimal dimension of work
call zgges(jobvsl,jobvsr,sort,selahattin,n,s,lda,t,ldb, &
             sdim, alpha, beta, vsl,ldvsl, vsr, ldvsr,work, &
			 lwork,rwork,bwork,info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))

! Calculate the sorted generalized Schur form. 
call zgges(jobvsl,jobvsr,sort,selahattin,n,s,lda,t,ldb, &
           sdim, alpha, beta, vsl,ldvsl, vsr, ldvsr,work, &
           lwork,rwork,bwork,info)

! Checking for the right number of stable eigenvalues

if ((abs(alpha(k))<abs(beta(k))) .or. (abs(alpha(k+1))>abs(beta(k+1)))) then
 print*, 'Error: wrong number of stable eigenvalues'
 retco = 1
 return
end if

z11  = vsr(1:k,1:k)
z21  = vsr(k+1:n,1:k)
s11i = s(1:k,1:k)
t11  = t(1:k,1:k)
call dlincg(k,z11,k,z11i,k)  ! Inverting z11
call dlincg(k,s11i,k,s11i,k) ! Inverting s11

f = real(matmul(z21,z11i))  ! f = z21*z11^(-1)
p = real(matmul(z11,s11i))  ! p = z11*s11^(-1)*t11*z11^(-1)
p = real(matmul(p,real(matmul(t11,z11i))))

end subroutine

! Selects stable generalized eigenvalues

logical function selahattin(sii,tii)
 complex*16 sii,tii
 selahattin = abs(sii)>abs(tii)
end function


end module