eigen Subroutine

  subroutine eigen (n, matrix, eigvec, eigval)
    !+The subroutine eigen  determines all eigenvalues and (if desired)
    !+all eigenvectors of a real square  n * n  matrix via the QR method
    !+in the version of Martin, Parlett, Peters, Reinsch and Wilkinson.
    !+
    !+###Literature
    !+1. Peters, Wilkinson: Eigenvectors of real and complex
    !+   matrices by LR and QR triangularisations,
    !+   Num. Math. 16, p.184-204, (1970); [PETE70]; contribution
    !+   II/15, p. 372 - 395 in [WILK71].
    !+2. Martin, Wilkinson: Similarity reductions of a general
    !+   matrix to Hessenberg form, Num. Math. 12, p. 349-368,(1968)
    !+   [MART 68]; contribution II,13, p. 339 - 358 in [WILK71].
    !+3. Parlett, Reinsch: Balancing a matrix for calculations of
    !+   eigenvalues and eigenvectors, Num. Math. 13, p. 293-304,
    !+   (1969); [PARL69]; contribution II/11, p.315 - 326 in
    !+   [WILK71].
    !+
    !+###Input parameters
    !+   n     :  integer; ( n > 0 )
    !+            size of matrix, number of eigenvalues
    !+
    !+   mat   :  n x n matrix;
    !+            input matrix
    !+
    !+###Output parameters
    !+   eivec :  n x n matrix;     (only if vec = 1)
    !+            matrix, if  vec = 1  that holds the eigenvectors
    !+            thus :
    !+            If the jth eigenvalue of the matrix is real then the
    !+            jth column is the corresponding real eigenvector;
    !+            if the jth eigenvalue is complex then the jth column
    !+            of eivec contains the real part of the eigenvector
    !+            while its imaginary part is in column j+1.
    !+            (the j+1st eigenvector is the complex conjugate
    !+            vector.)
    !+
    !+   valre :  vector of size n;
    !+            Real parts of the eigenvalues.
    !+
    !+   valim :  vector of size n;
    !+            Imaginary parts of the eigenvalues
    !+
    !+   cnt   :  Integer vector of size n;
    !+            vector containing the number of iterations for each
    !+            eigenvalue. (for a complex conjugate pair the second
    !+            entry is negative).
    integer      ,intent(in)                 :: n ! nlevels
    real(double) ,intent(in),dimension(n,n)  :: matrix
    real(double) ,intent(out),dimension(n,n) :: eigvec
    real(double) ,intent(out),dimension(n)   :: eigval
    real(double)     :: mat(0:n-1,0:n-1)
    real(double)     :: eivec(0:n-1,0:n-1)
    real(double)     :: valre(0:n-1) !real parts of eigenvalues
    real(double)     :: valim(0:n-1) !imaginary parts of eigenvalues
    integer          :: rc             !return code
    integer          :: cnt(0:n-1)   !Iteration counter
    integer          :: high, low
    real(double)     :: d(0:n-1), scale(0:n-1)
    integer          :: perm(0:n-1)

    cnt=0 ; d=0.d0
    mat(0:n-1,0:n-1)=matrix(1:n,1:n)
    !balance mat for nearly
    call balance(n, mat, scale, low, high)      !equal row and column
    !reduce mat to upper
    call elmhes(n, low, high, mat, perm)        !reduce mat to upper
    !Hessenberg form
    call elmtrans(n, low, high, mat, perm, eivec)
    !QR algorithm for eigenvalues and eigenvectors
    call hqr2(n, low, high, mat, valre, valim, eivec, cnt,rc)
    !reverse balancing to determine eigenvectors
    call balback(n, low, high, scale, eivec)
    if (rc.ne.0) then
      print*, 'matrix = '
      print*, matrix
      stop 'problem in eigen!'
    endif
    eigval(1:n)=valre(0:n-1)
    eigvec(1:n,1:n)=eivec(0:n-1,0:n-1)
  end subroutine eigen