C***********************************************************************
C
C                                   EIGCS
C
C PURPOSE:
C   Diagonalises a complex symmetric matrix, i.e. computes its
C   eigenvalues and right eigenvectors.
C DEPENDENCIES:
C   Needs LAPACK routine ZGEEV.
C NOTE:
C   An eigenvalue procedure for GENERAL complex matrices is used in this
C   routine. The performance will thus be rather poor if the matrices
C   are large.
C INPUT PARAMETERS:
C   n: dimension of "a", "d", and "u"; must be >= 1
C   a: complex symmetric matrix a to be diagonalised (full storage
C      mode, but only the upper triangular part is referenced)
C   lda: leading dimension of "a"; must be >= n
C   ldu: leading dimension of "u"; must be >= n
C OUTPUT PARAMETERS:
C   a: is overwritten
C   d: complex vector with the eigenvalues of "a"
C   u: complex matrix with the eigenvectors of "a" (stored columnwise)
C   err: error flag having the meaning
C        0: everything o.k.
C        1: invalid dimension of matrices
C        2: diagonalisation failed
C OTHER PARAMETERS:
C   caux: auxiliary complex array, must be >= 2*n
C   raux: auxiliary real array, must be >= 2*n
C
C MB 4/00
C
C***********************************************************************

      subroutine eigcs (n, a, lda, ldu, d, u, err, caux, raux)

      implicit none

      integer    n, lda, ldu, err
      real*8     raux(2*n)
      complex*16 a(lda,n), d(n), u(ldu,n), caux(2*n)

      integer    i, j
      complex*16 zero, one, cdummy(1), norm
      parameter (zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0))

C --- Initialise variables ---

      err = 0
      cdummy(1) = zero

C --- Check dimensions of arrays ---

      if (n.lt.1 .or. lda.lt.n .or. ldu.lt.n) then
         err = 1
         return
      endif

C --- Diagonalise matrix using LAPACK routine ZGEEV ---

      call zgeev('N','V',n,a,lda,d,cdummy,1,u,ldu,caux,2*n,raux,err)
      if (err .ne. 0) then
         err = 2
         return
      endif

C --- Re-normalise eigenvectors such that "u" is orthogonal ---

      do i = 1,n
         norm = zero
         do j = 1,n
            norm = norm+u(j,i)**2
         enddo
         norm = one/sqrt(norm)
         do j = 1,n
            u(j,i) = norm*u(j,i)
         enddo
      enddo

      return
      end
