C HAPACK_SHEIG.F
C   Gateway function for skew-Hamiltonian eigencomputations.
C
C Matlab call:  
C                    [...] = hapack_sheig(task,A,QG,...)
C
C   task = 1 :    [Ar,QGr] = hapack_sheig(1,A,QG)
C           [U1,U2,Ar,QGr] = hapack_sheig(1,A,QG)
C
C   task = 2 :    [Ar,QGr] = hapack_sheig(2,A,QG)
C           [U1,U2,Ar,QGr] = hapack_sheig(2,A,QG)
C
C   task = 3 : [scale,ilo] = hapack_sheig(3,A,QG,balanc)
C                 [Ar,QGr] = hapack_sheig(3,A,QG,balanc)
C       [Ar,QGr,scale,ilo] = hapack_sheig(3,A,QG,balanc)
C
C   task = 4 :    [Ar,QGr] = hapack_sheig(4,A,QG,select,lower)
C           [U1,U2,Ar,QGr] = hapack_sheig(4,A,QG,select,lower)
C           [U1,U2,Ar,QGr] = hapack_sheig(4,U1,U2,A,QG,select,lower)
C
C   task = 5 :     [wr,wi] = hapack_sheig(5,A,QG,balanc)
C       [Ar,QGr,wr,wi,sep] = hapack_sheig(5,A,QG,balanc)
C     [V1,V2,Ar,QGr,wr,wi] = hapack_sheig(5,A,QG,balanc)
C [V1,V2,Ar,QGr,wr,wi,sep] = hapack_sheig(5,A,QG,balanc)
C  [V1,V2,W1,W2,Ar,QGr,...
C                   wr,wi] = hapack_sheig(5,A,QG,balanc)
C  [V1,V2,W1,W2,Ar,QGr,...
C               wr,wi,sep] = hapack_sheig(5,A,QG,balanc)
C  [V1,V2,W1,W2,Ar,QGr,...
C           wr,we,sen,sep] = hapack_sheig(5,A,QG,balanc)
C
C   task = 6 :
C [W1,W2,Ar,QGr,wr,wi,sep] = hapack_sheig(6,A,QG,balanc)
C  [V1,V2,W1,W2,Ar,QGr,...
C               wr,wi,sen] = hapack_sheig(6,A,QG,balanc)
C
C Purpose:
C   To compute eigenvalues and invariant subspaces of a 2*N-by-2*N
C   skew-Hamiltoninan matrix W = [A, G; Q, A']. On input as well as on
C   output skew-Hamiltonian matrices are stored in compressed format in
c   an N-by-N array Ax and an N-by-(N+1) array QGx.
C
C   task =  1: Computes the PVL form Wr = [Ar, Qr; 0, Ar] of W with Ar
C              in upper Hessenberg form. Optionally, N-by-N matrices U1
C              and U2 are computed so that U = [U1, U2; -U2, U1]
C              satisfies U'*U=I and W = U*Wr*U'.
C
C   task =  2: Computes the skew-Hamiltonian Schur form Wr =
C              [Ar, Qr; 0, Ar] of W with Ar in real Schur form.
C              Optionally, N-by-N matrices U1 and U2 are computed so
C              that U = [U1, U2; -U2, U1] satisfies U'*U=I and
C              W = U*Wr*U'.
C
C   task =  3: Finds an accurate similarity transformation T such that
C              Wr = T\W*T has, as nearly as possible, approximately
C              equal row and column norms. T is a permutation of a
C              diagonal matrix and symplectic. T is stored in an
C              n-vector SCALE as described in DSHBAL.
C
C   task =  4: Reorders the eigenvalues selected by the arrays SELECT
C              and LOWER to the top left part of a matrix W in
C              skew-Hamiltonian Schur form. Optionally, the
C              corresponding orthogonal symplectic transformation
C              matrix U = [U1, U2; -U2, U1] is computed or a given
C              matrix U = [U1, U2; -U2, U1] is post-multiplied
C              by this matrix.
C
C   task = 5   or
C   task = 6:  Computes the j-th eigenvalue of W as
C              wr(j)+sqrt(-1)*wi(j), and optionally, its structured
C              condition number 1/sen(j), (real/complex part of) the
C              corresponding v-eigenvector [V1(j);V2(j)] with structured
C              condition number 1/sep(j) (only available for real
C              eigenvectors), (real/complex part of) the corresponding
C              w-eigenvector [W1(j);W2(j)] and the balanced, skew-
C              Hamiltonian Schur form Wr.
C
C Input parameters: 
C   task   - integer option to determine the computation to perform as
C            described above.
C   A      - real N-by-N matrix. If task = 4 A must be in real Schur
C            form.
C   QG     - real N-by-(N+1) matrix containing the skew-symmetric
C            matrices Q and G in compressed format. If task = 4 the part
C            containing the matrix Q must be zero.
C   U1,U2  - real N-by-N matrices.
C   balanc - determines whether W should be permuted (balanc = 'p'),
C            scaled (balanc = 's') or permuted and scaled (balanc = 'b')
C            prior to eigenvalue computations. Otherwise balanc = 'n'.
C   lower  - logical N-vector controlling which copy of each selected
C            eigenvalue is to be reordered to the top part of a
C            skew-Hamiltonian Schur form.
C   select - logical N-vector specifying the eigenvalues to be reordered
C            to the top part of a skew-Hamiltonian Schur form as
C            described in DSHSEN.
C
C Output parameters:
C   U1,U2  - real N-by-N matrices so that U = [U1, U2; -U2, U1]
C            satisfies U'*U=I and W = U*Wr*U'.
C   V1,V2  - real N-by-N matrices containing the v-eigenvectors of W.
C   W1,W2  - real N-by-N matrices containing the w-eigenvectors of W.
C   Ar     - real N-by-N matrix containing the (1,1) block of the
C            returned skew-Hamiltonian matrix Wr.
C   QGr    - real N-by-(N+1) matrix containing the skew-symmetric (2,1)
C            and (1,2) blocks of the returned skew-Hamiltonian matrix Wr
C            in compressed format.
C   ilo    - ilo-1 is the number of deflated eigenvalues in the balanced
C            skew-Hamiltonian matrix. 
C   scale  - real N-vector containing the scaling factors as returned by
C            DSHBAL.
C   sen    - real N-vector containing the reciprocals of the structured
C            eigenvalue condition numbers.
C   sep    - real N-vector containing the reciprocals of the structured
C            condition numbers for real v-eigenvectors.
C   wr     - real N-vector containing the real parts of the eigenvalues
C            of W.
C   wi     - real N-vector containing the imaginary parts of the
C            eigenvalues of W.
C
C Contributor:
C   D. Kressner, UCL CESAME, Louvain-la-Neuve, June 2003.
C
C **********************************************************************
C
      SUBROUTINE MEXFUNCTION( NLHS, PLHS, NRHS, PRHS )
      IMPLICIT NONE
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C
C     .. Mex-file interface parameters ..
      INTEGER           PLHS(*), PRHS(*)
      INTEGER           NLHS, NRHS
C
C     .. Mex-file integer functions ..
      INTEGER           mxCreateFull, mxGetPr
      INTEGER           mxGetM, mxGetN, mxGetString, mxIsChar, 
     $                  mxIsComplex, mxIsNumeric
C
C     .. Scalar parameters used by HAPACK subroutines ..
      CHARACTER         BALANC, JOBV, JOBW, SENSE
      INTEGER           INFO, LDA, LDQG, LDWORK, LDU1, LDU2, LDV1, LDV2,
     $                  LDW1, LDW2, N
      DOUBLE PRECISION  SBNRM
C
C     .. Allocatable arrays ..
C     !Fortran 90/95 (Fixed dimensions should be used with Fortran 77.)
      LOGICAL,          ALLOCATABLE :: LOWER(:), LWORK(:), SELECT(:)
      INTEGER,          ALLOCATABLE :: IWORK(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:), CS(:), DWORK(:), QG(:,:),
     $                                 SCALE(:), SEN(:), SEP(:), TAU(:),
     $                                 U1(:,:), U2(:,:), V1(:,:),
     $                                 V2(:,:), W1(:,:), W2(:,:), WI(:),
     $                                 WR(:)
C
C     .. Local variables and constant dimension arrays ..
      CHARACTER*120     TEXT
      LOGICAL           COMPSC, COMPSN, COMPSP, COMPU, COMPV, COMPW,
     $                  UPDU
      INTEGER           IERR, I, ILO, IP, M, NN, NB, NBMIN, NLM, TASK
      DOUBLE PRECISION  TEMP
C
C     .. External Functions ..
      INTEGER           ILAHAP
      EXTERNAL          ILAHAP
C
C     .. External subroutines ..
      EXTERNAL          DHAORD, DLACPY, DLASET, DOSGPV, DSHBAL, DSHES,
     $                  DSHEVX, DSHPVB
C
C     ..Intrinsic functions..
      INTRINSIC         DBLE, MAX, MIN
C
C     Check for proper number of arguments.
C
      IF ( NRHS.LT.3 ) THEN
         CALL mexErrMsgTxt
     $        ( 'HAPACK_SHEIG requires at least 3 input arguments' )
      ELSE IF ( NLHS.GT.10 ) THEN
         CALL mexErrMsgTxt
     $        ( 'HAPACK_SHEIG requires at most 10 output arguments' )
      END IF
C
C     Check dimensions of input parameters and read/set scalar
C     parameters.
C
      IF ( mxGetM( PRHS(1) ).NE.1 .OR. mxGetN( PRHS(1) ).NE.1 )
     $   CALL mexErrMsgTxt( 'TASK must be a scalar' )
      IF ( mxIsNumeric( PRHS(1) ).EQ.0 .OR. 
     $     mxIsComplex( PRHS(1) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'TASK must be an integer scalar' )
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(1) ), TEMP, 1 )
      TASK = TEMP
C
      IF ( TASK.LT.1 .OR. TASK.GT.6 )
     $   CALL mexErrMsgTxt
     $        ( 'The only admissible values of TASK are 1, 2, ..., 6')
C
      COMPSC = .TRUE.
      COMPSN = .FALSE.
      COMPSP = .FALSE.
      COMPU  = .FALSE.
      COMPV  = .FALSE.
      COMPW  = .FALSE.
      UPDU   = .FALSE.
      IF ( TASK.LE.2 ) THEN
         IF ( NLHS.GT.2 )
     $      COMPU = .TRUE.
         NLM = 4
      ELSE IF ( TASK.EQ.3 ) THEN
         IF ( NLHS.LT.2 )
     $      COMPSC = .FALSE.
         NLM = 4
      ELSE IF ( TASK.EQ.4 ) THEN
         IF ( NLHS.GT.2 )
     $      COMPU = .TRUE.
         IF ( NRHS.GE.7 )
     $      UPDU = .TRUE.
         NLM = 4
      ELSE IF ( TASK.EQ.5 ) THEN
         IF ( NLHS.LE.2 )
     $      COMPSC = .FAlSE.
         IF ( NLHS.GT.2 .AND. NLHS.NE.6 .AND. NLHS.NE.8 )
     $      COMPSP = .TRUE.
         IF ( NLHS.GT.5 ) 
     $      COMPV = .TRUE.
         IF ( NLHS.GE.8 )
     $      COMPW = .TRUE.
         IF ( NLHS.EQ.10 )
     $      COMPSN = .TRUE.
         NLM = 10
      ELSE IF ( TASK.EQ.6 ) THEN
         COMPW = .TRUE.
         IF ( NLHS.EQ.7 )
     $      COMPSP = .TRUE.
         IF ( NLHS.GT.7 ) THEN
            COMPV  = .TRUE.
            COMPSN = .TRUE.
         END IF
         NLM = 9
      END IF
C
      IF ( NLHS.GT.NLM ) THEN
         WRITE( TEXT, '('' HAPACK_SHEIG requires at most '',I4,
     $            '' output arguments'')' ) NLM
         CALL mexErrMsgTxt( TEXT )
      END IF
C
      IF ( TASK.LE.2 .AND. NRHS.GT.3 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at most 3 input arguments' )
      ELSE IF ( ( TASK.EQ.3.OR.TASK.GE.5 ) .AND. NRHS.GT.4 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at most 4 input arguments' )
      ELSE IF ( TASK.EQ.4 .AND. .NOT.COMPU .AND. NRHS.GT.5 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at most 5 input arguments' )
      ELSE IF ( TASK.EQ.4 .AND. NRHS.GT.7 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at most 7 input arguments' )
      END IF
C
      IF ( TASK.LE.2 .AND. NRHS.LT.3 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at least 3 input arguments' )
      ELSE IF ( ( TASK.EQ.3.OR.TASK.GE.5 ) .AND. NRHS.LT.4 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at least 4 input arguments' )
      ELSE IF ( TASK.EQ.4 .AND. NRHS.LT.5 ) THEN
         CALL mexErrMsgTxt
     $           ( 'HAPACK_SHEIG requires at least 5 input arguments' )
      END IF
C
C     Check BALANC.
C
      IP = 4
      IF ( TASK.EQ.3 .OR. TASK.GE.5 ) THEN
         IF ( mxGetM( PRHS(IP) ).NE.1 .OR. mxGetN( PRHS(IP) ).NE.1
     $        .OR. mxIsChar( PRHS(IP) ).NE.1 )
     $      CALL mexErrMsgTxt( 'BALANC must be a character' )
         IERR = mxGetString( PRHS(IP), BALANC, 1 )
      END IF
      IF ( UPDU ) THEN
         IP = 4
      ELSE
         IP = 2
      END IF
C
C     Check A.
C
      M = mxGetM( PRHS(IP) )
      N = mxGetN( PRHS(IP) )
      IF ( M.LT.0 .OR. N.LT.0 .OR. M.NE.N )
     $   CALL mexErrMsgTxt( 'A must be a square matrix' )
      IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $     mxIsComplex( PRHS(IP) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'A must be a real matrix' )
      IP = IP + 1
C
C     Check QG.
C
      M  = mxGetM( PRHS(IP) )
      NN = mxGetN( PRHS(IP) )
      IF ( M.NE.N .OR. NN.NE.N+1 )
     $   CALL mexErrMsgTxt( 'QG must be an n-by-n+1 matrix' )
      IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $     mxIsComplex( PRHS(IP) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'QG must be a real matrix' )
C
C     Check U1 and U2.
C
      IF ( UPDU ) THEN
         IP = 2
         M  = mxGetM( PRHS(IP) )
         NN = mxGetN( PRHS(IP) )
         IF ( M.NE.N .OR. NN.NE.N )
     $      CALL mexErrMsgTxt( 'U1 must be an n-by-n matrix' )
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $        mxIsComplex( PRHS(IP) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'U1 must be a real matrix' )
         IP = IP + 1
         M  = mxGetM( PRHS(IP) )
         NN = mxGetN( PRHS(IP) )
         IF ( M.NE.N .OR. NN.NE.N )
     $      CALL mexErrMsgTxt( 'U2 must be an n-by-n matrix' )
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $        mxIsComplex( PRHS(IP) ).EQ.1 )
     $      CALL mexErrMsgTxt( 'U2 must be a real matrix' )
         IP = 5
      END IF
C
C     Check SELECT AND LOWER.
C
      IF ( TASK.EQ.4 ) THEN
         IP = IP + 1
         M  = mxGetM( PRHS(IP) )
         NN = mxGetN( PRHS(IP) )
         IF ( .NOT.( ( M.EQ.1.AND.NN.EQ.N ).OR.
     $               ( M.EQ.N.AND.NN.EQ.1 ) ) )
     $      CALL mexErrMsgTxt( 'SELECT must be a vector of length n' )
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $        mxIsComplex( PRHS(IP) ).EQ.1 )
     $      CALL mexErrMsgTxt( 'SELECT must be a real vector' )
         IP = IP + 1
         M  = mxGetM( PRHS(IP) )
         NN = mxGetN( PRHS(IP) )
         IF ( .NOT.( ( M.EQ.1.AND.NN.EQ.N ).OR.
     $               ( M.EQ.N.AND.NN.EQ.1 ) ) )
     $      CALL mexErrMsgTxt( 'LOWER must be a vector of length n' )
         IF ( mxIsNumeric( PRHS(IP) ).EQ.0 .OR. 
     $        mxIsComplex( PRHS(IP) ).EQ.1 )
     $      CALL mexErrMsgTxt( 'LOWER must be a real vector' )
      END IF
C
C     Determine the lenghts of working arrays.
C
      LDA = MAX( 1, N )
      LDQG = MAX( 1, N )
      IF ( COMPU ) THEN
         LDU1 = MAX( 1, N )
         LDU2 = MAX( 1, N )
      ELSE
         LDU1 = 1
         LDU2 = 1
      END IF
      IF ( COMPV ) THEN
         LDV1 = MAX( 1, N )
         LDV2 = MAX( 1, N )
      ELSE
         LDV1 = 1
         LDV2 = 1
      END IF
      IF ( COMPW ) THEN
         LDW1 = MAX( 1, N )
         LDW2 = MAX( 1, N )
      ELSE
         LDW1 = 1
         LDW2 = 1
      END IF
C
C     Workspace requirements to enable blocked algorithms.
C
      LDWORK = 1
      IF ( TASK.EQ.1 ) THEN
         LDWORK = MAX( LDWORK, N )
         NB = ILAHAP( 1, 'DSHPVB', '', N, 1, -1 )
         NBMIN = MAX( 2, ILAHAP( 2, 'DSHPVB', '', N, 1, -1 ) )
         IF ( NB.GE.NBMIN )
     $      LDWORK = MAX( LDWORK, 8*N*NB + 3*NB )
         IF ( COMPU ) THEN
            LDWORK = MAX( LDWORK, 2*N )
            NB = ILAHAP( 1, 'DOSGSB', 'NN', N, N, N )
            NBMIN = MAX(  2, ILAHAP( 2, 'DOSGSB', 'NN', N, N, N ) )
            IF ( NB.GE.NBMIN )
     $         LDWORK = MAX( LDWORK, 8*N*NB + 15*NB*NB )
         END IF
      ELSE IF ( TASK.EQ.2 ) THEN
         NB = ILAHAP( 1, 'DSHPVB', '', N, 1, -1 )
         NBMIN = MAX( 2, ILAHAP( 2, 'DSHPVB', '', N, 1, -1 ) )
         IF ( NB.GE.NBMIN )
     $      LDWORK = MAX( LDWORK, 4*N + 8*N*NB + 3*NB )
         IF ( COMPU ) THEN
            LDWORK = MAX( LDWORK, (N+5)*N )
            NB = ILAHAP( 1, 'DOSMSB', 'NN', N, N, N )
            NBMIN = MAX(  2, ILAHAP( 2, 'DOSMSB', 'NN', N, N, N ) )
            IF ( NB.GE.NBMIN )
     $         LDWORK = MAX( LDWORK, (N+4)*N + 9*N*NB + 15*NB*NB )
         ELSE
            LDWORK = MAX( LDWORK, 5*N, (N+1)*N, 2*N*N-N )
         END IF
      ELSE IF ( TASK.EQ.4 ) THEN
         LDWORK = MAX( LDWORK, N )
      ELSE IF ( TASK.GE.5 ) THEN
         LDWORK = MAX( LDWORK, 5*N )
         NB = ILAHAP( 1, 'DSHPVB', '', N, 1, -1 )
         NBMIN = MAX( 2, ILAHAP( 2, 'DSHPVB', '', N, 1, -1 ) )
         IF ( NB.GE.NBMIN )
     $      LDWORK = MAX( LDWORK, 4*N + 8*N*NB + 3*NB )
         NB = ILAHAP( 1, 'DOSMSB', 'NN', N, N, N )
         NBMIN = MAX( 2, ILAHAP( 2, 'DOSMSB', 'NN', N, N, N ) )
         IF ( NB.LT.NBMIN )
     $      NB = 0
         IF ( COMPV.AND.COMPW ) THEN
            LDWORK = MAX( LDWORK, 5*N + 9*N*NB + 15*NB*NB )
         ELSE IF ( COMPV ) THEN
            LDWORK = MAX( LDWORK, N*N + 5*N + 9*N*NB + 15*NB*NB )
         ELSE IF ( COMPW ) THEN
            LDWORK = MAX( LDWORK, 2*N*N + 7*N + 9*N*NB + 15*NB*NB )
         END IF
         IF ( COMPSP ) THEN
            LDWORK = MAX( LDWORK, 3*N*N + 6*N )
         END IF
      END IF
C
C     Allocate variable dimension local arrays.
C     !Fortran 90/95
C
      ALLOCATE ( A(LDA,N), DWORK(LDWORK), QG(LDQG,N+1) )
      IF ( TASK.EQ.1 ) THEN
         ALLOCATE ( CS(2*N), TAU(N) )
      ELSE IF ( TASK.EQ.2 ) THEN
         ALLOCATE ( WI(N), WR(N) )
      ELSE IF ( TASK.EQ.3 ) THEN
         ALLOCATE ( SCALE(N) )
      ELSE IF ( TASK.EQ.4 ) THEN
         ALLOCATE ( LOWER(N), SELECT(N), WI(N), WR(N) )
      ELSE IF ( TASK.GE.5 ) THEN
         ALLOCATE ( IWORK(2*N), LWORK(2*N), SCALE(N), WI(N), WR(N) )
      END IF
      IF ( COMPU )
     $   ALLOCATE( U1(LDU1,N), U2(LDU2,N) )
      IF ( COMPV )
     $   ALLOCATE( V1(LDV1,N), V2(LDV2,N) )
      IF ( COMPW )
     $   ALLOCATE( W1(LDW1,N), W2(LDW2,N) )
      IF ( COMPSN )
     $   ALLOCATE( SEN(N) )
      IF ( COMPSP )
     $   ALLOCATE( SEP(N) )
C
C     Copy inputs from MATLAB workspace to locally allocated arrays.
C
      IP = 1
      IF ( UPDU ) THEN
         IP = IP + 1
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), U1, N*N )
         IP = IP + 1
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), U2, N*N )
      ELSE IF ( COMPU ) THEN
         CALL DLASET( 'All', N, N, ZERO, ONE, U1, LDU1 )
         CALL DLASET( 'All', N, N, ZERO, ZERO, U2, LDU2 )
      END IF
      IP = IP + 1
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), A, N*N )
      IP = IP + 1
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), QG, N*(N+1) )
      IF ( TASK.EQ.4 ) THEN
         IP = IP + 1
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), DWORK, N )
         DO 10  I = 1, N
             SELECT(I) = DWORK(I).EQ.ONE
   10    CONTINUE
         IP = IP + 1
         CALL mxCopyPtrToReal8( mxGetPr( PRHS(IP) ), DWORK, N )
         DO 20  I = 1, N
            LOWER(I) = DWORK(I).EQ.ONE
   20    CONTINUE
      END IF
C
C     Do the actual computations.
C
      IF ( TASK.EQ.1 ) THEN
         CALL DSHPVB( N, 1, A, LDA, QG, LDQG, CS, TAU, DWORK, LDWORK,
     $                INFO )
         IF ( INFO.EQ.0 ) THEN
            IF ( COMPU ) THEN
               CALL DLACPY( 'Lower', N, N, A, LDA, U1, LDU1 )
               CALL DLACPY( 'Lower', N, N, QG, LDQG, U2, LDU2 )
               CALL DOSGPV( N, 1, U1, LDU1, U2, LDU2, CS, TAU, DWORK,
     $                   LDWORK, INFO )
            END IF
            CALL DLASET( 'Lower', MAX(N-2,0), MAX(N-2,0), ZERO, ZERO,
     $                   A(3,1), LDA )
            CALL DLASET( 'Lower', MAX(N-1,0), MAX(N-1,0), ZERO, ZERO,
     $                   QG(2,1), LDQG )
         END IF
      ELSE IF ( TASK.EQ.2 ) THEN
         IF ( COMPU ) THEN
            CALL DSHES( 'U', N, A, LDA, QG, LDQG, U1, LDU1, U2, LDU2,
     $                  WR, WI, DWORK, LDWORK, INFO )
         ELSE
            CALL DSHES( 'N', N, A, LDA, QG, LDQG, DWORK, 1, DWORK, 1,
     $                  WR, WI, DWORK, LDWORK, INFO )
         END IF
      ELSE IF ( TASK.EQ.3 ) THEN
         CALL DSHBAL( BALANC, N, A, LDA, QG, LDQG, ILO, SCALE, INFO )
      ELSE IF ( TASK.EQ.4 ) THEN
         IF ( COMPU ) THEN
            CALL DHAORD( 'Skew-Hamiltonian', 'Update', SELECT, LOWER, N,
     $                   A, LDA, QG(1,2), LDQG, U1, LDU1, U2, LDU2, WR,
     $                   WI, M, DWORK, LDWORK, INFO )
         ELSE
            CALL DHAORD( 'Skew-Hamiltonian', 'No Update', SELECT, LOWER,
     $                   N, A, LDA, QG(1,2), LDQG, DWORK, 1, DWORK, 1,
     $                   WR, WI, M, DWORK, LDWORK, INFO )
         END IF
      ELSE IF ( TASK.GE.5 ) THEN
         IF ( COMPV ) THEN
            JOBV = 'V'
         ELSE
            JOBV = 'N'
         END IF
         IF ( COMPW ) THEN
            JOBW = 'W'
         ELSE
            JOBW = 'N'
         END IF
         IF ( COMPSN.AND.COMPSP ) THEN
            SENSE = 'B'
         ELSE IF ( COMPSN ) THEN
            SENSE = 'E'
         ELSE IF ( COMPSP ) THEN
            SENSE = 'V'
         ELSE
            SENSE = 'N'
         END IF
         CALL DSHEVX( BALANC, JOBV, JOBW, SENSE, N, A, LDA, QG,
     $                LDQG, WR, WI, V1, LDV1, V2, LDV2, W1, LDW1,
     $                W2, LDW2, ILO, SCALE, SBNRM, SEN, SEP,
     $                LWORK, IWORK, DWORK, LDWORK, INFO )
      END IF
C
C     Copy output to MATLAB workspace.
C
      IP = 1
      IF ( INFO.EQ.0 ) THEN
         IF ( COMPU ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( U1, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPU .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( U2, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPV ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( V1, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPV .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( V2, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPW .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( W1, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPW .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( W2, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPSC .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N, 0 )
            CALL mxCopyReal8ToPtr( A, mxGetPr( PLHS(IP) ), N*N )
            IP = IP + 1
         END IF
         IF ( COMPSC .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, N+1, 0 )
            CALL mxCopyReal8ToPtr( QG, mxGetPr( PLHS(IP) ), N*(N+1) )
            IP = IP + 1
         END IF
         IF ( TASK.EQ.3 .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, 1, 0 )
            CALL mxCopyReal8ToPtr( SCALE, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
         IF ( TASK.EQ.3 .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( 1, 1, 0 )
            TEMP = DBLE( ILO )
            CALL mxCopyReal8ToPtr( TEMP, mxGetPr( PLHS(IP) ), 1 )
            IP = IP + 1
         END IF
         IF ( TASK.GE.5 .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, 1, 0 )
            CALL mxCopyReal8ToPtr( WR, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
         IF ( TASK.GE.5 .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, 1, 0 )
            CALL mxCopyReal8ToPtr( WI, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
         IF ( COMPSN .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, 1, 0 )
            CALL mxCopyReal8ToPtr( SEN, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
         IF ( COMPSP .AND. NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( N, 1, 0 )
            CALL mxCopyReal8ToPtr( SEP, mxGetPr( PLHS(IP) ), N )
            IP = IP + 1
         END IF
      END IF
C
C     Deallocate local arrays.
C     !Fortran 90/95
C
      DEALLOCATE ( A, DWORK, QG )
      IF ( TASK.EQ.1 ) THEN
         DEALLOCATE ( CS, TAU )
      ELSE IF ( TASK.EQ.2 ) THEN
         DEALLOCATE ( WI, WR )
      ELSE IF ( TASK.EQ.3 ) THEN
         DEALLOCATE ( SCALE )
      ELSE IF ( TASK.EQ.4 ) THEN
         DEALLOCATE ( LOWER, SELECT, WI, WR )
      ELSE IF ( TASK.GE.5 ) THEN
         DEALLOCATE ( IWORK, LWORK, SCALE, WI, WR )
      END IF
      IF ( COMPU )
     $   DEALLOCATE( U1, U2 )
      IF ( COMPV )
     $   DEALLOCATE( V1, V2 )
      IF ( COMPW )
     $   DEALLOCATE( W1, W2 )
      IF ( COMPSN )
     $   DEALLOCATE( SEN )
      IF ( COMPSP )
     $   DEALLOCATE( SEP )
C
C     Error and warning handling.
C
      IF ( INFO.NE.0 ) THEN
         IF ( TASK.EQ.1 ) THEN
            IF ( COMPU ) THEN
               WRITE( TEXT, '('' INFO = '',I4,
     $                '' ON EXIT FROM DSHPVB/DOSGPV'')' )  INFO
            ELSE
               WRITE( TEXT,
     $                '('' INFO = '',I4,'' ON EXIT FROM DSHPVB'')' )
     $                INFO
            END IF
         ELSE IF ( TASK.EQ.2 ) THEN
            WRITE( TEXT, '('' INFO = '',I4,'' ON EXIT FROM DSHES'')' )
     $             INFO
         ELSE IF ( TASK.EQ.3 ) THEN
            WRITE( TEXT, '('' INFO = '',I4,'' ON EXIT FROM DSHBAL'')' )
     $             INFO
         ELSE IF ( TASK.EQ.4 ) THEN
            WRITE( TEXT, '('' INFO = '',I4,'' ON EXIT FROM DSHSEN'')' )
     $             INFO
         ELSE IF ( TASK.GE.5 ) THEN
            WRITE( TEXT, '('' INFO = '',I4,'' ON EXIT FROM DSHEVX'')' )
     $             INFO
         END IF
         CALL mexErrMsgTxt( TEXT )
      END IF
C
      RETURN
C *** Last line of HAPACK_SHEIG ***
      END