Purpose
 =======
 
      LA_GEES computes for a real/complex square matrix A, the 
 eigenvalues, the real-Schur/complex-Schur form T , and, optionally, the
 matrix of Schur vectors Z, where Z is orthogonal/unitary. This gives the
 Schur factorization
                        A = Z*T*Z^H.
 Optionally, it also orders the eigenvalues on the diagonal of the Schur 
 form so that selected eigenvalues are at the top left. The leading 
 columns of Z then form an orthonormal basis for the invariant subspace
 corresponding to the selected eigenvalues.
      A real matrix is in real-Schur form if it is block upper triangular
 with 1 by 1 and 2 by 2 blocks along the main diagonal. 2 by 2 blocks are
 standardized in the form
                       [ a  b ]
 	               [ c  a ]
 where b*c < 0. The eigenvalues of such a block are a +/- Sqrt(b*c).
 A complex matrix is in complex-Schur form if it is upper triangular.
 
 =========
 
        SUBROUTINE LA_GEES( A, <w>, VS=vs, SELECT=select, &
                                   SDIM=sdim, INFO=info )
           <type>(<wp>), INTENT(INOUT) :: A(:,:)
           <type>(<wp>), INTENT(OUT) :: <w(:)>
           <type>(<wp>), INTENT(OUT), OPTIONAL :: VS(:,:)
           INTERFACE
               LOGICAL FUNCTION SELECT(<w(j)>)
                  <type>(<wp>), INTENT(IN) :: <w(j)>
               END FUNCTION SELECT
           END INTERFACE
           OPTIONAL :: SELECT
           INTEGER, INTENT(OUT), OPTIONAL :: SDIM, INFO
        where
           <type> ::= REAL | COMPLEX
           <wp>   ::= KIND(1.0) | KIND(1.0D0)
           <w>    ::= WR, WI | W
           <w(:)> ::= WR(:), WI(:) | W(:)
           <w(j)> ::= WR(j) , WI(j) |  W(j)
 
 Arguments
 =========
 
 A        (input/output) REAL or COMPLEX square array, shape (:,:).
          On entry, the matrix A.
          On exit, the Schur form T.
 <w>      (output) REAL or COMPLEX array, shape (:) with size(w) = 
          size(A,1).
          The computed eigenvalues in the order in which they appear on 
          the diagonal of the Schur form T.
          <w(:)> ::= WR(:), WI(:) | W(:),
          where
          WR(:), WI(:) are of REAL type (for the real and imaginary 
 	  parts) and W(:) is of COMPLEX type.
          Note: If A is real, then a complex-conjugate pair appear 
          consecutively, with the eigenvalue having the positive
 	  imaginary part appearing first.
 VS       Optional (output) REAL or COMPLEX square array, shape (:,:) 
          with size(VS,1) = size(A,1).
          The matrix Z of Schur vectors.
 SELECT   Optional (input) LOGICAL FUNCTION.
          LOGICAL FUNCTION SELECT( <w(j)> )
            <type>(<wp>), INTENT(IN) :: <w(j)>
          where
            <type> ::= REAL | COMPLEX
            <wp>   ::= KIND(1.0) | KIND(1.0D0)
            <w(j)> ::= WR(j) , WI(j) | W(j)
          1. SELECT must be declared as EXTERNAL or as an explicit 
             interface in the calling (sub)program.
          2. SELECT is called by LA_GEES for every computed eigenvalue
 	    w(j) (but only once for a complex conjugate pair when A is
 	    real). It is used to select the eigenvalues that will be
            ordered to the top left of the Schur form. The eigenvalue
 	    w(j) is selected if SELECT(w(j)) has the value .TRUE.
          3. A selected complex eigenvalue may no longer satisfy 
 	    SELECT(w(j)) = .TRUE. after ordering, since ordering may 
 	    change the value of complex eigenvalues (especially if the
 	    eigenvalue is ill-conditioned). In this case INFO is set to
 	    size(A,1) + 2 (see INFO below).
          Note: Select must be present if SDIM is desired.
 SDIM     Optional (output) INTEGER.
          The number of eigenvalues (after sorting) for which 
 	   SELECT=.TRUE. (If A is real, complex conjugate pairs for which 
 	   SELECT=.TRUE. for either eigenvalue count as 2).
 INFO     Optional (output) INTEGER.
          = 0: successful exit.
          < 0: if INFO = -i, the i-th argument had an illegal value.
          > 0: if INFO = i, and i is 
              <= n: the QR algorithm failed to compute all the 
                    eigenvalues; elements 1:ilo-1 and i+1:n of w contain
                    those eigenvalues which have converged. VS contains 
                    the matrix which reduces A to its partially converged
                    Schur form.
              = n+1: the eigenvalues could not be reordered because some 
                    eigenvalues were not sufficiently separated (the 
 		   problem is very ill-conditioned).
              = n+2: after reordering, some leading complex eigenvalues 
 	           in the Schur form no longer satisfy SELECT = .TRUE. 
                   This can be caused by ordinary roundoff or underflow
                   due to scaling.
              n is the order of A.
          If INFO is not present and an error occurs, then the program is
          terminated with an error message.