      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
