      SUBROUTINE MVZERO (L,M,N,NMAX,LMAX,MAX,A,B,C,D,G,H,NZER,NRANK,
     +                   ZR,ZI,BETA,SUM,DUMMY)
C
C     ***** PURPOSE:
C      THIS SUBROUTINE COMPUTES THE FINITE MULTIVARIABLE ZEROS OF THE
C      LINEAR SYSTEM
C
C               X'(T) = A*X(T) + B*U(T)
C                Y(T) = C*X(T) + D*U(T)
C
C      WHERE A IS N X N, B IS N X M, C IS L X N, AND D IS L X M. THE
C      ALGORITHM EMPLOYED IS DESCRIBED BELOW UNDER ALGORITHM NOTES.
C
C     ***** CALLING SEQUENCE PARAMETERS:
      INTEGER L,M,N,NMAX,LMAX,MAX,NZER,NRANK
      DOUBLE PRECISION A(NMAX,N),B(NMAX,M),C(LMAX,N),D(LMAX,M),G(MAX,1),
     +                 H(MAX,1),ZR(1),ZI(1),BETA(1),SUM(1),DUMMY(1)
C
C     ***** LOCAL VARIABLES:
      INTEGER I,IERR,J
      DOUBLE PRECISION EPS,EPSP1,EPS1,GNORM,T
      LOGICAL MATZ
C
C     ***** FORTRAN FUNCTIONS CALLED:
      DOUBLE PRECISION DABS
C
C     ***** SUBROUTINES CALLED:
C     ZEROS,QZHES(EISPACK),QZIT(EISPACK),QZVAL(EISPACK)
C
C     ------------------------------------------------------------------
C
C     ***** PARAMETER DESCRIPTION:
C     ON INPUT:
C
C        L,M,N      INTEGERS
C                   NUMBERS OF OUTPUTS, INPUTS, AND STATES, RESPEC-
C                   TIVELY, OF THE LINEAR SYSTEM WHOSE ZEROS ARE TO BE
C                   DETERMINED.
C
C        NMAX,LMAX  INTEGERS
C                   ROW DIMENSIONS OF THE ARRAYS CONTAINING A,B AND C,D,
C                   RESPECTIVELY, AS DECLARED IN THE MAIN CALLING
C                   PROGRAM DIMENSION STATEMENT.
C
C        MAX        INTEGER
C                   ROW DIMENSION OF THE (N+L) X (N+M) ARRAYS G,H AS
C                   DECLARED IN THE MAIN CALLING PROGRAM DIMENSION
C                   STATEMENT.
C                   *** NOTE THAT MAX MUST BE  .GE. (N+L)
C
C        A,B,C,D    DOUBLE PRECISION (NMAX,N),(NMAX,M),(LMAX,N),(LMAX,M)
C                   ARRAYS CONTAINING THE SYSTEM MATRICES; NOTE THAT
C                   THESE ARRAYS ARE NOT MODIFIED BY THE SUBROUTINE.
C
C        G,H        DOUBLE PRECISION (MAX,1),(MAX,1)
C                   SCRATCH ARRAYS OF SIZE (N+L) X (N+M).
C
C        BETA,SUM,  DOUBLE PRECISION (1),(1),(1)
C          DUMMY    SCRATCH VECTORS OF LENGTHS N, MAX(L,M), MAX(L,M,N),
C                   RESPECTIVELY.  NOTE: LENGTH N WILL SUFFICE FOR ALL
C                   THREE VECTORS ASSUMING THERE ARE AT LEAST AS MANY
C                   STATES AS INPUTS OR OUTPUTS.
C
C     ON OUTPUT:
C
C        NZER       INTEGER
C                   = K   ON NORMAL RETURN WITH K (FROM 0 TO N)
C                         FINITE MULTIVARIABLE ZEROS.
C                   = -K  ON ABNORMAL RETURN IF THE QZ ALGORITHM
C                         APPLIED TO THE FINITE PENCIL FAILED TO
C                         CONVERGE AT ROW K.
C
C        NRANK      INTEGER
C                   THE NORMAL RANK OF THE SYSTEM TRANSFER FUNCTION;
C                   GENERICALLY, THIS NUMBER WILL BE MIN(L,M);  IF
C                   NRANK IS RETURNED AS 0 THE SYSTEM IS SAID TO BE
C                   DEGENERATE BUT MAY STILL HAVE MEANINGFUL ZEROS (CF.
C                   REFERENCES (1) AND (2)).
C
C        ZR,ZI      DOUBLE PRECISION (1),(1)
C                   VECTORS (OF LENGTH AT MOST N) CONTAINING THE REAL
C                   AND IMAGINARY PARTS, RESPECTIVELY, OF THE NZER
C                   FINITE MULTIVARIABLE ZEROS.
C
C     ***** ALGORITHM NOTES:
C      1. THE ZEROS ARE DETERMINED BY MEANS OF THE QZ ALGORITHM
C         (SEE REF. (3)) APPLIED TO A FINITE MATRIX PENCIL WHICH
C         IS DETERMINED AS A RESULT OF AN ORTHOGONAL STRUCTURE
C         ALGORITHM APPLIED TO THE ORIGINAL SYSTEM PENCIL;  REFS.
C         (1),(2), AND (3) SHOULD ALL BE CONSULTED FOR FURTHER
C         ALGORITHMIC DETAILS.
C      2. TWO IMPORTANT ARRAYS ARE PRODUCED BY THE QZ ALGORITHM:
C         THE ALPHA'S WHICH ARE COMPLEX AND THE BETA'S WHICH ARE
C         REAL.  THE ZEROS ARE DETERMINED AS RATIOS (ALPHA/BETA)
C         OF THESE NUMBERS.  THE BETA'S ARE AVAILABLE IF DESIRED
C         IN THE SCRATCH VECTOR BETA AND THE ALPHA'S MAY BE
C         RECOVERED BY MULTIPLYING COMPONENTS OF ZR AND ZI BY THE
C         RESPECTIVE COMPONENTS OF BETA.
C
C     ***** REFERENCES:
C      (1) EMAMI-NAEINI, A., AND P. VANDOOREN, COMPUTATION OF ZEROS
C           OF LINEAR MULTIVARIABLE SYSTEMS, REPORT NO. NA-80-03,
C           DEP'T. OF COMP. SCI., STANFORD UNIVERSITY, JULY, 1980.
C
C      (2) LAUB, A.J., AND B.C. MOORE, CALCULATION OF TRANSMISSION
C           ZEROS USING QZ TECHNIQUES, AUTOMATICA, 14(1978),557-566.
C
C      (3) MOLER, C.B., AND G.W. STEWART, AN ALGORITHM FOR GENER-
C           ALIZED MATRIX EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL.,
C           10(1973),241-256.
C
C     ***** SOFTWARE HISTORY:
C      WRITTEN BY    ALAN J. LAUB
C                    DEP'T. OF ELEC. ENGRG. - SYSTEMS
C                    UNIV. OF SOUTHERN CALIFORNIA
C                    LOS ANGELES, CA 90089
C                    PH.: (213) 743-5535
C      JANUARY, 1981.
C      MOST RECENT VERSION:  AUG. 3, 1983.
C
C     ------------------------------------------------------------------
C
      MATZ = .FALSE.
      EPS1 = -1.0D0
C
C     DETERMINE MACHINE PRECISION
C
      EPS = 1.0D0
   10 CONTINUE
      EPS = EPS/2.0D0
      EPSP1 = EPS+1.0D0
      IF (EPSP1 .GT. 1.0D0) GO TO 10
      EPS = 2.0D0*EPS
C
C     COMPUTE   SQRT(EPS) * 1-NORM OF (A B)
C                                     (C D)
C
      GNORM = 0.0D0
      DO 40 J = 1,N
         T = 0.0D0
         DO 20 I = 1,N
            T = T+DABS(A(I,J))
   20    CONTINUE
         DO 30 I = 1,L
            T = T+DABS(C(I,J))
   30    CONTINUE
         IF (T .GT. GNORM) GNORM = T
   40 CONTINUE
      DO 70 J = 1,M
         T = 0.0D0
         DO 50 I = 1,N
            T = T+DABS(B(I,J))
   50    CONTINUE
         DO 60 I = 1,L
            T = T+DABS(D(I,J))
   60    CONTINUE
         IF (T .GT. GNORM) GNORM = T
   70 CONTINUE
      EPS = DSQRT(EPS)*GNORM
C
C     GET REDUCED ORDER PENCIL WHICH WILL CONTAIN THE FINITE ZEROS
C
      CALL ZEROS (A,B,C,D,M,NMAX,N,LMAX,L,MAX,EPS,H,G,
     +            NZER,NRANK,SUM,DUMMY)
C
C     APPLY QZ ALGORITHM TO THE FINITE PENCIL
C
      IF (NZER .EQ. 0) RETURN
      IF (NZER .EQ. 1) GO TO 100
      CALL QZHES (MAX,NZER,G,H,MATZ,G)
      CALL QZIT (MAX,NZER,G,H,EPS1,MATZ,G,IERR)
      IF (IERR .EQ. 0) GO TO 80
      NZER = -IERR
      RETURN
   80 CONTINUE
      CALL QZVAL (MAX,NZER,G,H,ZR,ZI,BETA,MATZ,G)
C
C     COMPUTE THE MULTIVARIABLE ZEROS
C
      DO 90 I = 1,NZER
         ZR(I) = ZR(I)/BETA(I)
         ZI(I) = ZI(I)/BETA(I)
   90 CONTINUE
      RETURN
  100 CONTINUE
      ZR(1) = G(1,1)/H(1,1)
      ZI(1) = 0.0D0
      RETURN
      END
