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