SUBROUTINE ORDER (A,B,Z,NMAX,N,EPS,IFAIL,TYPE,IFIRST,IND)
C
C     *****PARAMETERS:
      INTEGER IFAIL,IFIRST,N,NMAX,IND(N)
      DOUBLE PRECISION EPS,A(NMAX,N),B(NMAX,N),Z(NMAX,N)
      LOGICAL TYPE
C
C     *****LOCAL VARIABLES:
      INTEGER I,II,III,IS,ISTEP,K,L,LS,LS1,LS2,L1,L2,NUM
C     DOUBLE PRECISION EPS2
C
C     *****FORTRAN FUNCTIONS:
      DOUBLE PRECISION DABS
C
C     *****SUBROUTINES CALLED:
C     EXCHQZ
C
C     ---------------------------------------------------------------
C
C     *****PURPOSE:
C     GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG
C     MATRIX A WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS SUBROUTINE
C     REORDERS THE DIAGONAL BLOCKS ALONG WITH THE GENERALIZED
C     EIGENVALUES CORRESPONDING TO THE REGULAR MATRIX PENCIL
C     A - LAMBDA*B  BY CONSTRUCTING ROW AND COLUMN EQUIVALENCE
C     TRANSFORMATIONS QT AND ZT.  THE COLUMN TRANSFORMATIONS ARE
C     THEN APPLIED TO THE MATRIX Z.
C
C     REF.:  VAN DOOREN, P., A GENERALIZED EIGENVALUE APPROACH FOR
C     SOLVING RICCATI EQUATIONS, SIAM J. SCI. STAT. COMPUT.,
C     VOL. 2, NO. 2, JUNE 1981, 121-135.
C
C     *****PARAMETER DESCRIPTION:
C
C     ON INPUT:
C
C       NMAX    INTEGER
C               ROW DIMENSION OF THE ARRAYS CONTAINING A,B,Z AS
C               DECLARED IN THE MAIN CALLING PROGRAM DIMENSION STATEMENT;
C
C       N       INTEGER
C               ORDER OF THE MATRICES A,B,Z;
C
C       A       REAL(NMAX,N)
C               UPPER HESSENBERG MATRIX WITH 1X1 OR 2X2 DIAGONAL
C               BLOCKS.  ELEMENTS OUTSIDE THE UPPER HESSENBERG
C               STRUCTURE ARE ARBITRARY;
C
C       B       REAL(NMAX,N)
C               UPPER TRIANGULAR MATRIX.  ELEMENTS OUTSIDE THE
C               UPPER TRIANGULAR STRUCTURE ARE ARBITRARY;
C
C       EPS     REAL
C               REQUIRED ABSOLUTE ACCURACY OF THE RESULTS.  NORMALLY
C               EQUAL TO THE MACHINE PRECISION;
C
C       TYPE    LOGICAL
C               CONTROL PARAMETER THAT SPECIFIES THE REGIONS OF THE
C               THE COMPLEX PLANE THAT THE GENERALIZED EIGENVALUES
C               ARE ORDERED BY.  TO CONTROL THE REGION THAT APPEARS
C               FIRST, SEE IFIRST BELOW
C               = .TRUE.  GENERALIZED EIGENVALUES ARE ORDERED BY
C                    REGION INSIDE THE COMPLEX LEFT HALF PLANE OR
C                    OUTSIDE THIS REGION
C               = .FALSE.  GENERALIZED EIGENVALUES ORDERED BY REGION
C                    INSIDE THE UNIT CIRCLE OR OUTSIDE THIS REGION;
C
C       IFIRST  INTEGER
C               CONTROL PARAMETER THAT SPECIFIES WHICH OF THE REGIONS
C               SPECIFIED BY TYPE(SEE ABOVE) APPEARS FIRST(I.E. IN THE
C               UPPER LEFT NXN BLOCK)
C               = -1  INSIDE REGION APPEARS FIRST
C               = +1  OUTSIDE REGION APPEARS FIRST
C               IFIRST IS ALTERED BY THIS SUBROUTINE;
C
C       IND     INTEGER(N)
C               WORKING ARRAY THAT IS ALTERED BY THIS SUBROUTINE.
C
C     ON OUTPUT:
C
C       A,B     UPPER HESSENBERG MATRIX, UPPER TRIANGULAR MATRIX
C               REORDERED AS SPECIFIED BY TYPE AND IFIRST(SEE ABOVE);
C
C       Z       REAL(NMAX,N)
C               THIS ARRAY IS OVERWRITTEN BY THE PRODUCT OF THE
C               CONTENTS OF THE ARRAY Z(UPON ENTRY INTO THIS
C               SUBROUTINE), AND THE COLUMN TRANSFORMATIONS ZT
C               (CALCULATED BY THIS SUBROUTINE);
C
C       IFAIL   INTEGER
C               ERROR FLAG
C               = 1  INDICATES ATTEMPTED REORDERING FAILED
C               = 0  NORMAL RETURN.
C
C     *****ALGORITHM NOTES:
C     NONE
C
C     *****HISTORY:
C     ORIGINAL VERSION THAT SORTED BY UNIT CIRCLE REGION OF COMPLEX
C     PLANE WRITTEN BY P. VAN DOOREN("A GENERALIZED EIGENVALUE
C     APPROACH FOR SOLVING RICCATI EQUATIONS", INTERNAL REPORT
C     NA-80-02, DEPT. OF COMPUTER SCIENCE, STANFORD UNIVERSITY,
C     1980).  THIS VERSION MODIFIED BY W. F. ARNOLD(DEPT. OF
C     ELECTRICAL ENGINEERING - SYSTEMS, UNIV. OF SOUTHERN CALIF.,
C     LOS ANGELES, CA  90089) TO INCLUDE THE SORTING CONTROL PARAMETER
C     "TYPE", SEPT 1982.
C
C     ----------------------------------------------------------------
C
      IFAIL = 1
      NUM = 0
      L = 0
      LS = 1
C
C--JDB/CASCADE: PATCH TO MAKE ALGORITHM LESS SENSITIVE TO INEXACT ZEROS, 7/2/87
C
C     EPS2 = SQRT(EPS)
C
C     DETERMINE SIZE AND REGION LOCATION OF BLOCKS
C
   10 CONTINUE
      L = L+LS
      IF (L .GT. N) GO TO 50
      IS = -IFIRST
      L1 = L+1
      IF (L1 .GT. N) GO TO 20
C--JDB/CASCADE: REPLACED FOLLOWING LINE 6/16/87
      IF (A(L1,L) .EQ. 0.0D0) GO TO 20
C     IF (A(L1,L) .LE. EPS2*DMAX1(DABS(A(L,L)),DABS(A(L1,L1)))) GO TO 20
C
C     2X2 BLOCK
C
      LS = 2
      IF (.NOT. TYPE) GO TO 15
      IF ((A(L,L)*B(L1,L1)+A(L1,L1)*B(L,L)-A(L1,L)*B(L,L1))/
     X    (B(L,L)*B(L1,L1)) .LT. 0.0D0) IS = IFIRST
      GO TO 30
   15 CONTINUE
      IF (DABS(A(L,L)*A(L1,L1)-A(L1,L)*A(L,L1)) .LT.
     +    DABS(B(L,L)*B(L1,L1))) IS= IFIRST
      GO TO 30
C
C     1X1 BLOCK
C
   20 CONTINUE
      LS = 1
      IF (.NOT. TYPE) GO TO 25
      IF (A(L,L)*B(L,L) .LT. 0.0D0) IS = IFIRST
      IF (A(L,L) .LT. 0.0D0 .AND. B(L,L) .EQ. 0.0D0) IS = IFIRST
      GO TO 30
   25 CONTINUE
      IF (DABS(A(L,L)) .LT. DABS(B(L,L))) IS = IFIRST
   30 CONTINUE
      NUM = NUM+1
      IND(NUM) = LS*IS
      GO TO 10
C
C     REORDER BLOCKS
C
   50 CONTINUE
      L2 = 1
      I = 0
   60 CONTINUE
      I = I+1
      IF (IND(I) .GT. 0) GO TO 70
      L2 = L2-IND(I)
      GO TO 60
   70 CONTINUE
      K = I
   80 CONTINUE
      L2 = L2+IND(K)
   85 CONTINUE
      K = K+1
      IF (K .GT. NUM) GO TO 100
      IF (IND(K) .GT. 0) GO TO 80
C
C     INTERCHANGE BLOCK K BEFORE BLOCK I
C
      ISTEP = K-I
      LS2 = -IND(K)
      L = L2
      DO 90 II = 1,ISTEP
         III = K-II
         LS1 = IND(III)
         L = L-LS1
         CALL EXCHQZ (A,B,Z,NMAX,N,L,LS1,LS2,EPS,IFAIL)
         IF (IFAIL .EQ. 1) RETURN
         IND(III+1) = IND(III)
   90 CONTINUE
      IND(I) = -LS2
      I = I+1
      L2 = L2+LS2
      GO TO 85
  100 CONTINUE
      IFAIL = 0
      RETURN
C
C     LAST LINE OF ORDER
C
      END