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