SUBROUTINE DFRMG (FIRST,NA,NB,NC,L,M,N,A,B,C,FREQ,AINVB,G,H,IPVT, + ZC,ZR,EVRE,EVIM,IERR,INFO,RCOND) INTEGER NA,NB,NC,L,M,N,IPVT(N),IERR,INFO DOUBLE PRECISION A(NA,N),B(NB,M),C(NC,N),ZR(N),EVRE(N),EVIM(N), + RCOND DOUBLE COMPLEX FREQ,AINVB(NB,M),G(NC,M),H(NA,N),ZC(N) LOGICAL FIRST C C *** PURPOSE: C C DFRMG TAKES REAL MATRICES A (N X N), B (N X M), AND C (L X N) C AND FORMS THE COMPLEX FREQUENCY RESPONSE MATRIX C G(FREQ) := C * (((FREQ * I) - A)-INVERSE) * B C WHERE I = (N X N) IDENTITY MATRIX AND FREQ IS A COMPLEX C SCALAR PARAMETER TAKING VALUES ALONG THE IMAGINARY AXIS FOR C CONTINUOUS-TIME SYSTEMS AND ON THE UNIT CIRCLE FOR DISCRETE- C TIME SYSTEMS. C C ON ENTRY: C C FIRST LOGICAL C SET = .TRUE. FOR THE FIRST CALL OF DFRMG WHEREUPON C IT IS SET TO .FALSE. FOR ALL SUBSEQUENT CALLS; MAY BE C SET = .FALSE. INITIALLY IF A IS ALREADY IN UPPER C HESSENBERG FORM AND NEITHER BALANCING NOR EIGENVALUES C OF A ARE DESIRED. C C NA INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY A C (AND THE COMPLEX ARRAY H) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C NB INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY B C (AND THE COMPLEX ARRAY AINVB) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C NC INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY C C (AND THE COMPLEX ARRAY G) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C L INTEGER C THE NUMBER OF ROWS OF C (THE NUMBER OF OUTPUTS). C C M INTEGER C THE NUMBER OF COLUMNS OF B (THE NUMBER OF INPUTS). C C N INTEGER C THE ORDER OF THE MATRIX A (THE NUMBER OF STATES); C ALSO = NUMBER OF COLUMNS OF C = NUMBER OF ROWS OF B. C C A DOUBLE PRECISION(NA,N) C A REAL N X N MATRIX (THE SYSTEM MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C B DOUBLE PRECISION(NB,M) C A REAL N X M MATRIX (THE INPUT MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C C DOUBLE PRECISION(NC,N) C A REAL L X N MATRIX (THE OUTPUT MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C FREQ DOUBLE COMPLEX C A COMPLEX SCALAR (THE FREQUENCY PARAMETER). C C IERR INTEGER C --FIRST = .FALSE. : -IERR IS IGNORED. C --FIRST = .TRUE. : -IF IERR .GE. 0, A WILL BE BALANCED C AND THE EIGENVALUES OF A WILL BE C COMPUTED (RECOMMENDED OPTION). C -IF IERR .LT. 0, A WILL NOT BE BAL- C ANCED AND THE EIGENVALUES OF A WILL C NOT BE COMPUTED. C C INFO INTEGER C --IF INFO .LT. 0, THE CONDITION OF (FREQ*I - A) WITH C RESPECT TO INVERSION WILL BE ESTIMATED (RECOMMENDED C OPTION). C --IF INFO .GE. 0, THE CONDITION OF (FREQ*I - A) WITH C RESPECT TO INVERSION WILL NOT BE ESTIMATED. C C ON RETURN: C C G DOUBLE COMPLEX(NC,M) C THE FREQUENCY RESPONSE MATRIX G(FREQ). C C A,B,C A IS IN UPPER HESSENBERG FORM WHILE B AND C HAVE BEEN C SUITABLY TRANSFORMED; IF FIRST .EQ. .FALSE. THESE C ARRAYS ARE NOT FURTHER MODIFIED. C C EVRE, DOUBLE PRECISION(N),DOUBLE PRECISION(N) C EVIM THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE C EIGENVALUES OF THE MATRIX A IF DFRMG WAS INITIALLY C CALLED WITH FIRST = .TRUE. .AND. IERR .GE. 0. C C IERR IF (FIRST = .TRUE. .AND. IERR .GE. 0 ON INPUT) THEN C IERR IS OUTPUT AS THE IERR PARAMETER OF THE EISPACK C SUBROUTINE HQR (HQR DOCUMENTATION MAY BE CONSULTED FOR C DETAILS); NORMAL RETURN IS 0. C C INFO --IGNORE IF INFO .LT. 0 ON INPUT. C --IF INFO .GE. 0 ON INPUT, THIS IS OUTPUT AS THE INFO C PARAMETER OF SUBROUTINE ZHEFA (ZHEFA MAY BE CONSULTED C FOR DETAILS); NORMAL RETURN IS 0. C C RCOND REAL C --IGNORE IF INFO .GE. 0 ON INPUT. C --IF INFO .LT. 0 ON INPUT, THIS IS OUTPUT AS THE RCOND C PARAMETER OF SUBROUTINE ZHECO (ZHECO MAY BE CONSULTED C FOR DETAILS); NORMAL RETURN IS WHEN C (1.0 + RCOND) .GT. 1.0. C C AINVB DOUBLE COMPLEX(NB,M) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C H DOUBLE COMPLEX(NA,N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C IPVT INTEGER(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C ZC DOUBLE COMPLEX(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C ZR DOUBLE PRECISION(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C THIS VERSION DATED JUNE 1982. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C BALANC(EISPACK),ZHECO,ZHEFA,ZHESL,HQR(EISPACK),DHETR C C INTERNAL VARIABLES: C INTEGER I,IGH,J,K,KK,KP,LOW DOUBLE PRECISION T DOUBLE COMPLEX ZDUM C C FORTRAN FUNCTIONS CALLED: C DOUBLE COMPLEX DCMPLX C C INTERNAL FUNCTION: C DOUBLE PRECISION DREAL DREAL(ZDUM) = ZDUM C IF (.NOT. FIRST) GO TO 150 IF (IERR .GE. 0) GO TO 10 LOW = 1 IGH = N GO TO 90 10 CONTINUE CALL BALANC (NA,N,A,LOW,IGH,ZR) C C ADJUST B AND C MATRICES BASED ON INFORMATION IN THE VECTOR C ZR WHICH DESCRIBES THE BALANCING OF A AND IS DEFINED IN THE C SUBROUTINE BALANC C DO 40 K = 1,N KK = N-K+1 IF (KK .GE. LOW .AND. KK .LE. IGH) GO TO 40 IF (KK .LT. LOW) KK = LOW-KK KP = ZR(KK) IF (KP .EQ. KK) GO TO 40 C C PERMUTE ROWS OF B C DO 20 J = 1,M T = B(KK,J) B(KK,J) = B(KP,J) B(KP,J) = T 20 CONTINUE C C PERMUTE COLUMNS OF C C DO 30 I = 1,L T = C(I,KK) C(I,KK) = C(I,KP) C(I,KP) = T 30 CONTINUE C 40 CONTINUE IF (IGH .EQ. LOW) GO TO 80 DO 70 K = LOW,IGH T = ZR(K) C C SCALE COLUMNS OF PERMUTED C C DO 50 I = 1,L C(I,K) = C(I,K)*T 50 CONTINUE C C SCALE ROWS OF PERMUTED B C DO 60 J = 1,M B(K,J) = B(K,J)/T 60 CONTINUE C 70 CONTINUE 80 CONTINUE 90 CONTINUE C C REDUCE A TO HESSENBERG FORM BY ORTHOGONAL SIMILARITIES AND C ACCUMULATE THE ORTHOGONAL TRANSFORMATIONS INTO B AND C C CALL DHETR (NA,NB,NC,L,M,N,LOW,IGH,A,B,C,ZR) IF (IERR .LT. 0) GO TO 140 C C TEMPORARILY STORE HESSENBERG FORM OF A IN THE ARRAY H C DO 110 J = 1,N DO 100 I = 1,N H(I,J) = DCMPLX(A(I,J),0.0D0) 100 CONTINUE 110 CONTINUE C C COMPUTE THE EIGENVALUES OF A IF THAT OPTION IS REQUESTED C CALL HQR (NA,N,LOW,IGH,A,EVRE,EVIM,IERR) C C RESTORE UPPER HESSENBERG FORM OF A TO THE ARRAY A C DO 130 J = 1,N DO 120 I = 1,N A(I,J) = DREAL(H(I,J)) 120 CONTINUE 130 CONTINUE 140 CONTINUE FIRST = .FALSE. IF (IERR .GT. 0) RETURN C C UPDATE H := (FREQ * I) - A WITH APPROPRIATE VALUE OF FREQ C 150 CONTINUE DO 170 J = 1,N DO 160 I = 1,N H(I,J) = -DCMPLX(A(I,J),0.0D0) 160 CONTINUE H(J,J) = FREQ+H(J,J) 170 CONTINUE IF (INFO .LT. 0) GO TO 180 C C FACTOR THE COMPLEX HESSENBERG MATRIX WITH NO CONDITION C ESTIMATE C CALL ZHEFA (H,NA,N,IPVT,INFO) IF (INFO .EQ. 0) GO TO 190 RETURN 180 CONTINUE C C FACTOR THE COMPLEX HESSENBERG MATRIX AND ESTIMATE ITS C CONDITION C CALL ZHECO (H,NA,N,IPVT,RCOND,ZC) T = 1.0D0+RCOND IF (T .EQ. 1.0D0) RETURN 190 CONTINUE C C COMPUTE (H-INVERSE)*B C DO 220 J = 1,M DO 200 I = 1,N ZC(I) = DCMPLX(B(I,J),0.0D0) 200 CONTINUE CALL ZHESL (H,NA,N,IPVT,ZC) DO 210 I = 1,N AINVB(I,J) = ZC(I) 210 CONTINUE 220 CONTINUE C C COMPUTE C*(H-INVERSE)*B C DO 260 J = 1,M DO 230 I = 1,L G(I,J) = (0.0D0,0.0D0) 230 CONTINUE DO 250 K = 1,N DO 240 I = 1,L G(I,J) = G(I,J)+DCMPLX(C(I,K),0.0D0)*AINVB(K,J) 240 CONTINUE C C G NOW CONTAINS THE DESIRED FREQUENCY RESPONSE MATRIX C 250 CONTINUE 260 CONTINUE RETURN END