c================================================================================
c  Diagonalize complex eigen equation AV = EBV
c================================================================================
      PROGRAM DiagonalizeCAB
      implicit real (a-h,o-z)
      implicit integer (i-n)
      parameter(nDIM=200)

      complex*16 A(nDIM,nDIM), B(nDIM,nDIM), V(nDIM,nDIM)
      real*8     E(nDIM), W(nDIM,7)
      integer*4  LW(nDIM), wd(nDIM,7)

      integer*4 nOrb, ier
      real*8    eps
      character*80 InFile, OutFile, CheckFile
      common /params/nOrb

      InFile    = 'matrix.in'
      OutFile   = 'matrix.out'
      CheckFile = 'matrix.txt'

      open(15, file=CheckFile, status='unknown')

      Call RdMtx(InFile, A, B, nOrb, eps);
      Write(15,*) 'nOrb=', nOrb, '   eps=', eps
      do 10 i = 1, nOrb
        do 20 j = 1, nOrb
          write(15, 100) i, j, A(i,j), B(i,j)
  100     format(i3,3x,i3,5x,4(g12.6,x))
   20   continue
   10 continue

      Call EigCH(A, nOrb, nDIM, nOrb, nOrb,
     1            epseig, W, LW, E, V, ier)
c      Call EigCAB(A, B, nOrb, nDIM, nOrb, nOrb,
c     1            epseig, W, E, V, ier, lw, wd)
      Write(15,*) 'Eigen values'
      do 110 i = 1, nOrb
          write(15, *) 'E=', E(i)
  110 continue
      do 220 i = 1, nOrb
        do 230 j = 1, nOrb
          write(15, 200) i, j, V(i,j)
  200     format(i3,3x,i3,5x,2(g12.6,x))
  230   continue
  220 continue

      Call SaveEig(OutFile, E, V, nOrb)

      close(15)

      END

c================================================================================
c
c================================================================================
      Subroutine SaveEig(OutFile, E, V, nOrb)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      parameter(nDIM=200)

      character*(*) OutFile
      real*8        E(*)
      complex*16    V(nDIM,nDIM)
      integer*4     nOrb

      integer i, j
      integer funit

      funit = 8
c      Write(*, *) 'Save to ', OutFile
      open(funit, file=OutFile, status='unknown')

      Write(funit,*) nOrb
      Write(funit,*) 'Eigen values'
      do 10 i = 1, nOrb
          write(funit, *) i, E(i)
   10 continue

      Write(funit,*) 'Eigen vectors'
      do 110 i = 1, nOrb
        do 120 j = 1, nOrb
          write(funit, *) i, j, real(V(i,j)), imag(V(i,j))
  120   continue
  110 continue

      close(funit)

      return
      End

c================================================================================
c
c================================================================================
      Subroutine RdMtx(InFile, A, B, nOrb, eps)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      parameter(nDIM=200)
      
      character*(*) InFile
      integer       funit
      complex*16,intent(out) :: A(nDIM,nDIM), B(nDIM,nDIM)
      integer*4, intent(out) :: nOrb
      real*8,    intent(out) :: eps

      character*80 s
      integer i, j, ii, jj
      real*8  re, img

      funit = 4
c      Write(*, *) 'Read ', InFile
      open(funit, file=InFile, status='old')

      read(funit, *) nOrb, eps
      read(funit, *) s
      do 10 i = 1, nOrb
        do 20 j = 1, nOrb
          read(funit, *) ii, jj, re, img
          if(ii .ne. i .or. jj .ne. j) then
            write(*, *) 'Error: Inconsistent (i,j): ', i, j, ii, jj, 
     -                  re, img
            stop
          endif
          A(i,j) = dcmplx(re, img)
   20   continue
   10 continue

      read(funit, *) s
      do 110 i = 1, nOrb
        do 120 j = 1, nOrb
          read(funit, *) ii, jj, re, img
          if(ii .ne. i .or. jj .ne. j) then
            write(*, *) 'Error: Inconsistent (i,j): ', i, j, ii, jj, 
     -                  re, img
            stop
          endif
          B(i,j) = dcmplx(re, img)
  120   continue
  110 continue

      close(funit)

      return
      End