      SUBROUTINE RESID(NR,NRR,NRW,NRX,N,M,E,A,B,CQC,R,S,RI,RSDM,X,W1,
     X                 W2,WK,IPVT,RTOL,EFLAG,RFLAG,RDFLG,SFLAG,RSD,
     X                 TYPE,NOUT)
C
C     *****PARAMETERS:
      INTEGER NR,NRR,NRW,NRX,N,M,IPVT(N),NOUT
      CHARACTER EFLAG,RFLAG,RDFLG,SFLAG
      DOUBLE PRECISION E(NR,N),A(NR,N),B(NR,M),CQC(NR,N),R(NR,M),
     X    S(NR,M),RI(NR,M),RSDM(NRR,N),X(NRX,N),W1(NRW,N),W2(NRW,N),
     X    WK(N),RTOL,RSD
      LOGICAL TYPE
C
C     *****LOCAL VARIABLES:
      INTEGER I
      DOUBLE PRECISION COND
C
C     *****FORTRAN FUNCTIONS:
C     NONE
C
C     *****FUNCTION SUBPROGRAMS:
      DOUBLE PRECISION D1NRM
C
C     *****SUBROUTINES CALLED:
C     MADD, MLINEQ, MMUL, MQF, MQFA, MSUB, MULA, TRNATB
C
C     --------------------------------------------------------------
C
C     *****PURPOSE:
C     THIS SUBROUTINE CALCULATES THE RESIDUAL MATRIX AND ITS 1-NORM
C     FOR THE GARE AS FOLLOWS:
C
C     CONTINUOUS:
C
C     RSDM = AT*X*E + ET*X*A - (BT*X*E+ST)T*RI*(BT*X*E+ST) + CQC
C
C
C     DISCRETE:
C
C     RSDM = AT*X*A - ET*X*E + CQC
C            - (BT*X*A+ST)T*((R+BT*X*B)**-1)*(BT*X*A+ST)
C
C
C     WHERE T DENOTES MATRIX TRANSPOSE.
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 ARRAYS CONTAINING THE MATRICES
C               E, A, B, CQC, R, RI AND S AS DECLARED IN THE MAIN
C               CALLING PROGRAM DIMENSION STATEMENT;
C
C       NRR     INTEGER
C               ROW DIMENSION OF THE ARRAY CONTAINING THE MATRIX
C               RSDM AS DECLARED IN THE MAIN CALLING PROGRAM
C               DIMENSION STATEMENT;
C
C       NRW     INTEGER
C               ROW DIMENSION OF THE ARRAYS CONTAINING THE MATRICES
C               W1 AND W2 AS DECLARED IN THE MAIN CALLING PROGRAM
C               DIMENSION STATEMENT;
C
C       NRX     INTEGER
C               ROW DIMENSION OF THE ARRAY CONTAINING THE MATRIX X
C               AS DECLARED IN THE MAIN CALLING PROGRAM DIMENSION
C               STATEMENT;
C
C       N       INTEGER
C               ORDER OF THE SQUARE MATRICES E, A, CQC, RSDM AND X
C               ROW DIMENSION OF THE MATRICES B AND S;
C
C       M       INTEGER
C               ORDER OF THE SQUARE MATRICES R AND RI
C               COLUMN DIMENSION OF THE MATRICES B AND S;
C
C       E       REAL(NR,N)
C               MODEL DESCRIPTOR MATRIX;
C
C       A       REAL(NR,N)
C               MODEL SYSTEM MATRIX;
C
C       B       REAL(NR,M)
C               MODEL INPUT MATRIX;
C
C       CQC     REAL(NR,N)
C               MATRIX PRODUCT CT*Q*C WHERE T DENOTES MATRIX
C               TRANSPOSE;
C
C       R       REAL(NR,M)
C               CONTROL WEIGHTING MATRIX;
C
C       S       REAL(NR,M)
C               STATE - INPUT CROSS-WEIGHTING MATRIX;
C
C       RI      REAL(NR,M)
C               INVERSE OF THE CONTROL WEIGHTING MATRIX;
C
C       X       REAL(NRX,N)
C               SOLUTION MATRIX FOR THE GARE WHOSE RESIDUAL IS TO BE
C               DETERMINED;
C
C       W1,W2   REAL(NRW,N)
C               SCRATCH ARRAYS OF SIZE AT LEAST N BY N;
C
C       WK      REAL(N)
C               WORK VECTOR OF LENGTH AT LEAST N;
C
C       IPVT    INTEGER(N)
C               WORK VECTOR OF LENGTH AT LEAST N;
C
C       RTOL    REAL
C               TOLERANCE ON THE CONDITION ESTIMATE OF R+BT*X*B WITH
C               RESPECT TO INVERSION (DISCRETE PROBLEM).  IF THIS
C               TOLERANCE IS EXCEEDED AN ERROR MESSAGE IS PRINTED
C               AND AN ERROR RETURN IS MADE;
C
C       EFLAG   CHARACTER
C               FLAG SET TO 'Y' IF E IS OTHER THAN THE IDENTITY
C               MATRIX;
C
C       RFLAG   CHARACTER
C               FLAG SET TO 'Y' IF R IS OTHER THAN THE IDENTITY
C               MATRIX;
C
C       RDFLG   CHARACTER
C               FLAG SET TO 'Y' IF R IS A DIAGONAL MATRIX;
C
C       SFLAG   CHARACTER
C               FLAG SET TO 'Y' IF S IS OTHER THAN THE ZERO MATRIX;
C
C       TYPE    LOGICAL
C               = .TRUE.  FOR CONTINUOUS-TIME SYSTEM
C               = .FALSE. FOR DISCRETE-TIME SYSTEM;
C
C       NOUT    INTEGER
C               UNIT NUMBER OF OUTPUT DEVICE FOR ERROR WARNING
C               MESSAGES.
C
C     ON OUTPUT:
C
C       RSDM    REAL(NRR,N)
C               THE RESIDUAL MATRIX CALCULATED AS INDICATED ABOVE;
C
C       IPVT(1) ERROR FLAG AS FOLLOWS
C               = 0  NORMAL RETURN
C               = 1  THE CONDITION ESTIMATE OF R+BT*X*B WITH RESPECT
C                    TO INVERSION EXCEEDS THE TOLERANCE RTOL;
C
C       RSD     REAL
C               THE 1-NORM OF THE RESIDUAL MATRIX DIVIDED BY THE
C               1-NORM OF THE SOLUTION MATRIX.
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
      CALL TRNATB(NR,NRR,N,M,B,RSDM)
      IF (TYPE) GO TO 80
C
C     DISCRETE-TIME CASE
C
      CALL MULA(NRR,NRX,M,N,N,RSDM,X,WK)
      CALL MMUL(NRR,NR,NRW,M,M,N,RSDM,B,W1)
      IF (RFLAG .EQ. 'Y' .OR. RFLAG .EQ. 'y') GO TO 20
      DO 10 I=1,M
         W1(I,I) = W1(I,I) + 1.0D0
   10 CONTINUE
      GO TO 50
   20 CONTINUE
      IF (RDFLG .EQ. 'Y' .OR. RDFLG .EQ. 'y') GO TO 30
      CALL MADD(NRW,NR,NRW,M,M,W1,R,W1)
      GO TO 50
   30 CONTINUE
      DO 40 I=1,M
         W1(I,I) = W1(I,I) + R(I,I)
   40 CONTINUE
   50 CONTINUE
      CALL MULA(NRR,NR,M,N,N,RSDM,A,WK)
      IF (SFLAG .NE. 'Y' .AND. SFLAG .NE. 'y') GO TO 60
      CALL TRNATB(NR,NRW,N,M,S,W2)
      CALL MADD(NRR,NRW,NRR,M,N,RSDM,W2,RSDM)
   60 CONTINUE
      CALL TRNATB(NRR,NRW,M,N,RSDM,W2)
      CALL MLINEQ(NRW,NRR,M,N,W1,RSDM,COND,IPVT,WK)
      IF (COND .GT. RTOL) GO TO 130
      CALL MULA(NRW,NRR,N,M,N,W2,RSDM,WK)
      CALL MQF(NRX,NR,NRW,N,N,X,A,W1,WK)
      IF (EFLAG .EQ. 'Y' .OR. EFLAG .EQ. 'y') GO TO 70
      CALL MSUB(NRW,NRX,NRW,N,N,W1,X,W1)
      GO TO 120
   70 CONTINUE
      CALL MQF(NRX,NR,NRR,N,N,X,E,RSDM,WK)
      CALL MSUB(NRW,NRR,NRW,N,N,W1,RSDM,W1)
      GO TO 120
   80 CONTINUE
C
C     CONTINUOUS-TIME CASE
C
      CALL TRNATB(NR,NRW,N,N,A,W2)
      IF (EFLAG .NE. 'Y' .AND. EFLAG .NE. 'y') GO TO 90
      CALL MMUL(NRX,NR,NRW,N,N,N,X,E,W1)
      CALL MULA(NRW,NRW,N,N,N,W2,W1,WK)
      CALL MULA(NRR,NRW,M,N,N,RSDM,W1,WK)
      GO TO 100
   90 CONTINUE
      CALL MULA(NRW,NRX,N,N,N,W2,X,WK)
      CALL MULA(NRR,NRX,M,N,N,RSDM,X,WK)
  100 CONTINUE
      IF (SFLAG .NE. 'Y' .AND. SFLAG .NE. 'y') GO TO 110
      CALL TRNATB(NR,NRW,N,M,S,W1)
      CALL MADD(NRR,NRW,NRR,M,N,RSDM,W1,RSDM)
  110 CONTINUE
      CALL TRNATB(NRW,NRW,N,N,W2,W1)
      CALL MADD(NRW,NRW,NRW,N,N,W2,W1,W1)
      IF (RFLAG .EQ. 'Y' .OR. RFLAG .EQ. 'y')
     X                 CALL MQF(NR,NRR,NRW,M,N,RI,RSDM,W2,WK)
      IF (RFLAG .NE. 'Y' .AND. RFLAG .NE. 'y')
     X                 CALL MQFA(NRR,NRW,M,N,RSDM,W2,WK)
  120 CONTINUE
      CALL MSUB(NRW,NRW,NRR,N,N,W1,W2,RSDM)
      CALL MADD(NRR,NR,NRR,N,N,RSDM,CQC,RSDM)
      RSD = D1NRM(NRR,N,N,RSDM) / D1NRM(NRX,N,N,X)
      RETURN
  130 CONTINUE
      WRITE(NOUT,140)
  140 FORMAT(//' R + BT*X*B  ILL-CONDITIONED WRT INVERSION'//)
      RETURN
C
C     LAST LINE OF RESID
C
      END
