      SUBROUTINE RICSOL(NR,NRD,NN,N,G,F,E,Z,ALFR,ALFI,BETA,CPERM,
     X                  CSCALE,IND,IORD,IBAL,TYPE,EFLAG)
C
C     *****PARAMETERS:
      INTEGER NR,NRD,NN,N,IND(NN),IORD,IBAL
      CHARACTER EFLAG
      DOUBLE PRECISION G(NRD,NN),F(NRD,NN),E(NR,N),Z(NRD,NN),
     X          ALFR(NN),ALFI(NN),BETA(NN),CPERM(NN),CSCALE(NN)
      LOGICAL TYPE
C
C     *****LOCAL VARIABLES:
      INTEGER I,IERR,IFAIL,IGH,J,LOW,NPJ
      DOUBLE PRECISION COND,EPS,EPS1
      LOGICAL MATZ
C
C     *****FORTRAN FUNCTIONS:
C     NONE
C
C     *****SUBROUTINES CALLED:
C     BALGBK, BALGEN, MLINEQ, MULB, ORDER, QZHESW, QZITW, QZVAL
C
C     --------------------------------------------------------------
C
C     *****PURPOSE:
C     GIVEN THE 2N BY 2N MATRIX PENCIL
C
C              LAMBDA * F  -  G
C
C     THIS SUBROUTINE FINDS THE ORTHOGONAL MATRIX Z SUCH THAT
C
C              Q * (LAMBDA * F  -  G) * Z
C
C     IS IN GENERALIZED ORDERED REAL SCHUR FORM.  FURTHERMORE, THE
C     UPPER LEFT N BY N BLOCK OF THE TRANSFORMED PENCIL CONTAINS
C     THE EIGENVALUES SPECIFIED BY THE PARAMETER IORD.  THE
C     SUBROUTINE THEN SOLVES THE LINEAR SYSTEM
C
C              X * E * Z11  =  Z21
C
C     FOR X, WHERE Z11 AND Z21 ARE THE UPPER AND LOWER LEFT N BY N
C     BLOCKS OF Z.
C
C     REF.:  ARNOLD, W.F., "ON THE NUMERICAL SOLUTION OF
C     ALGEBRAIC MATRIX RICCATI EQUATIONS," PHD THESIS, USC,
C     DECEMBER 1983.
C
C     *****PARAMETER DESCRIPTION:
C
C     ON INPUT:
C
C       NR      INTEGER
C               ROW DIMENSION OF THE ARRAY CONTAINING THE MATRIX E
C               AS DECLARED IN THE MAIN CALLING PROGRAM DIMENSION
C               STATEMENT;
C
C       NRD     INTEGER
C               ROW DIMENSION OF THE ARRAYS CONTAINING THE MATRICES
C               G, F AND Z AS DECLARED IN THE MAIN CALLING PROGRAM
C               DIMENSION STATEMENT;
C
C       NN      INTEGER
C               ORDER OF THE SQUARE MATRICES G AND F;
C
C       N       INTEGER
C               ORDER OF THE SQUARE MATRIX E;
C
C       G       REAL(NRD,NN)
C               PENCIL MATRIX CORRESPONDING TO THE GENERALIZED RICCATI
C               PROBLEM, WHOSE CONTENTS ARE ALTERED BY THIS ROUTINE;
C
C       F       REAL(NRD,NN)
C               PENCIL MATRIX CORRESPONDING TO THE GENERALIZED RICCATI
C               PROBLEM, WHOSE CONTENTS ARE ALTERED BY THIS ROUTINE;
C
C       E       REAL(NR,N)
C               DESCRIPTOR MATRIX
C
C       CPERM   REAL(NN)
C               WORKING VECTOR OF SIZE AT LEAST NN;
C
C       CSCALE  REAL(NN)
C               WORKING VECTOR OF SIZE AT LEAST NN;
C
C       IND     INTEGER(NN)
C               WORKING VECTOR OF SIZE AT LEAST NN;
C
C       IORD    INTEGER
C               PARAMETER SPECIFYING THE SPECTRUM OF THE UPPER LEFT N
C               BY N BLOCK OF THE ORDERED REAL SCHUR FORM AS FOLLOWS:
C               =  1  GENERALIZED EIGENVALUES WHOSE MAGNITUDE IS LESS
C                     THAN UNITY
C               =  0  ANY ORDER
C               = -1  GENERALIZED EIGENVALUES WHOSE REAL PARTS ARE
C                     LESS THAN ZERO;
C
C       IBAL    INTEGER
C               PARAMETER SPECIFYING THE BALANCING BEING EMPLOYED AS
C               FOLLOWS:
C               = 1  WARD BALANCING
C               = 0  NO BALANCING
C
C       TYPE    LOGICAL
C               = .TRUE.  FOR CONTINUOUS-TIME SYSTEM
C               = .FALSE. FOR DISCRETE-TIME SYSTEM;
C
C       EFLAG   CHARACTER
C               FLAG SET TO 'Y' IF E IS A DESCRIPTOR MATRIX THAT IS
C               OTHER THAN THE IDENTITY MATRIX.
C
C     ON OUTPUT:
C
C       F       CONTAINS THE SOLUTION MATRIX X COMPUTED AS SHOWN
C               ABOVE;
C
C       Z       REAL(NRD,NN)
C               CONTAINS THE MATRIX PRODUCT
C                    ( E  0 )     ( Z11  Z12)
C                    (      )  *  (         )
C                    ( 0  I )     ( Z21  Z22)
C               WHERE Z IS THE ORTHOGONAL TRANSFORMATION MATRIX
C               DESCRIBED ABOVE;
C
C       ALFR    REAL(NN)
C               REAL PARTS OF THE DIAGONAL ELEMENTS THAT WOULD RESULT
C               IF THE Q AND Z TRANSFORMATIONS WERE APPLIED TO THE
C               G MATRIX SUCH THAT IT WOULD BE REDUCED COMPLETELY TO
C               TRIANGULAR FORM AND THE DIAGONAL ELEMENTS OF THE
C               TRANSFORMED F MATRIX (ALSO TRIANGULAR) WOULD BE REAL
C               AND POSITIVE;
C
C       ALFI    REAL(NN)
C               IMAGINARY PARTS OF THE DIAGONAL ELEMENTS THAT WOULD
C               RESULT IF THE Q AND THE Z TRANSFORMATIONS WERE
C               APPLIED TO THE G MATRIX SUCH THAT IT WOULD BE REDUCED
C               COMPLETELY TO TRIANGULAR FORM AND THE DIAGONAL
C               ELEMENTS OF THE TRANSFORMED F MATRIX (ALSO TRIANGULAR)
C               WOULD BE REAL AND POSITIVE.  NONZERO VALUES OCCUR IN
C               PAIRS; THE FIRST MEMBER IS POSITIVE AND THE SECOND
C               MEMBER IS NEGATIVE;
C
C       BETA    REAL(NN)
C               REAL NONNEGATIVE DIAGONAL ELEMENTS OF F THAT WOULD
C               RESULT IF G WERE REDUCED COMPLETELY TO TRIANGULAR
C               FORM; THE GENERALIZED EIGENVALUES ARE THEN GIVEN BY
C               THE RATIOS ((ALFR + I*ALFI)/BETA);
C
C       CPERM(1)
C               CONDITION ESTIMATE OF E*Z11 WITH RESPECT TO INVERSION;
C
C       IND(1)  ERROR FLAG AS FOLLOWS
C               = 0  INDICATES NORMAL RETURN
C               = NONZERO IF MORE THAT 50 ITERATIONS WERE REQUIRED TO
C               DETERMINE THE DIAGONAL BLOCKS FOR THE QUASITRIANGULAR
C               FORM;
C
C       IND(2)  ERROR FLAG AS FOLLOWS
C               = 0  INDICATES NORMAL RETURN
C               = 1  INDICATES ATTEMPTED REORDERING FAILED.
C
C     *****ALGORITHM NOTES:
C     NONE
C
C     *****HISTORY:
C     THIS SUBROUTINE WAS WRITTEN BY W.F. ARNOLD, NAVAL WEAPONS
C     CENTER, CODE 35104, CHINA LAKE, CA  93555, AS PART OF THE
C     SOFTWARE PACKAGE RICPACK, SEPTEMBER 1983.
C     MODIFIED BY ALAN J. LAUB (UCSB): 01/06/85
C
C     --------------------------------------------------------------
C
      MATZ = .TRUE.
      LOW = 1
      IGH = NN
C
C     CALCULATE MACHINE PRECISION (EPS)
C
      EPS = 1.0D0
   10 CONTINUE
      EPS = EPS/2.0D0
      EPS1 = EPS + 1.0D0
      IF (EPS1 .GT. 1.0D0) GO TO 10
      EPS = EPS*2.0D0
C
C     WARD BALANCING IF REQUESTED
C
      IF (IBAL .EQ. 0) GO TO 20
      CALL BALGEN(NN,NRD,G,NRD,F,LOW,IGH,CSCALE,CPERM,Z)
   20 CONTINUE
C
C     CALCULATE GENERALIZED EIGENVALUES
C
      CALL QZHESW(NRD,NN,G,F,LOW,IGH,MATZ,Z)
      CALL QZITW(NRD,NN,G,F,LOW,IGH,EPS,MATZ,Z,IERR)
      CALL QZVAL(NRD,NN,G,F,ALFR,ALFI,BETA,MATZ,Z)
      EPS1 = F(NN,1)
C
C     ORDER EIGENVALUES
C
      IF(IORD .EQ. 0) GO TO 25
      CALL ORDER(G,F,Z,NRD,NN,EPS1,IFAIL,TYPE,IORD,IND)
   25 CONTINUE
      IF (IBAL .EQ. 0) GO TO 30
      CALL BALGBK(NN,NRD,Z,NN,LOW,IGH,CSCALE,CPERM)
   30 CONTINUE
C
C     CALCULATE RICCATI SOLUTION
C
      IF(EFLAG .NE. 'Y' .AND. EFLAG .NE. 'y') GO TO 35
      CALL MULB(NR,NRD,N,N,N,E,Z,CPERM)
   35 CONTINUE
      DO 50 J=1,N
         NPJ = N + J
         DO 40 I=1,N
            G(I,J) = Z(J,I)
            F(I,J) = Z(NPJ,I)
   40    CONTINUE
   50 CONTINUE
      CALL MLINEQ(NRD,NRD,N,N,G,F,COND,IND,CPERM)
      CPERM(1) = COND
      IND(1) = IERR
      IND(2) = IFAIL
      RETURN
C
C     LAST LINE OF RICSOL
C
      END
