SUBROUTINE ZEROS (A,B,C,D,M,NMAX,N,PMAX,P,MAX,EPS,BF,AF,NU,
     +                  RANK,SUM,DUMMY)
C
C        THIS ROUTINE EXTRACTS FROM THE SYSTEM MATRIX OF A STATE-
C        SPACE SYSTEM A:NXN, B:NXM, C:PXN, D:PXM  A REGULAR PENCIL
C        LAMBDA*BF - AF  WHERE BF AND AF ARE NU X NU WITH THE NU
C        INVARAIANT ZEROS OF THE SYSTEM AS GENERALIZED EIGENVALUES.
C        THE SUBROUTINE ZEROS REQUIRES THE SUBROUTINES RDUCE, HOUSH,
C        PIVOT, TR1, AND TR2.  THE PARAMETERS IN THE CALLING SEQUENCE
C        ARE:
C          INPUT:
C          A,B,C,D    THE SYSTEM DESCRIPTOR MATRICES
C          M,N,P      THE NUMBER OF INPUTS, STATES, AND OUTPUTS
C          PMAX,NMAX  THE FIRST DIMENSION OF C,D AND A,B RESPECTIVELY
C          MAX        THE FIRST DIMENSION OF AF,BF
C          EPS        THE ABSOLUTE TOLERANCE OF THE DATA (NOISE LEVEL);
C                     IT SHOULD BE LARGER THAN THE MACHINE PRECISION
C                     TIMES SOME NORM OF (A,B,C,D)
C          OUTPUT:
C          BF,AF      THE COEFFICIENT MATRICES OF THE REDUCED PENCIL
C          NU         THE NUMBER OF (FINITE) INVARIANT ZEROS
C          RANK       THE NORMAL RANK OF THE TRANSFER FUNCTION
C          WORKING SPACE:
C          SUM        A VECTOR OF DIMENSION AT LEAST MAX(M,P)
C          DUMMY      A VECTOR OF DIMENSION AT LEAST MAX(M,N,P)
C
C     *** PARTIALLY MODIFIED BY ALAN J. LAUB, UNIVERSITY OF SOUTHERN
C         CALIFORNIA, FEB. 1981 AND AUG. 1983.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL ZERO
      INTEGER P,PMAX,PP,RANK,RO,SIGMA
      DIMENSION A(NMAX,N),B(NMAX,M),C(PMAX,N),D(PMAX,M),AF(MAX,1),
     +          BF(MAX,1),SUM(1),DUMMY(1)
      MM = M
      NN = N
      PP = P
C
C     CONSTRUCT THE COMPOUND MATRIX  (B A) OF DIMEN. (N+P)X(M+N)
C                                    (D C)
C
      IF (MM .EQ. 0) GO TO 15
      DO 10 I = 1,NN
         DO 10 J = 1,MM
   10       BF(I,J) = B(I,J)
   15 DO 20 I = 1,NN
         DO 20 J = 1,NN
   20       BF(I,J+MM) = A(I,J)
      IF (PP .EQ. 0) GO TO 45
      IF (MM .EQ. 0) GO TO 35
      DO 30 I = 1,PP
         DO 30 J = 1,MM
   30       BF(I+NN,J) = D(I,J)
   35 DO 40 I = 1,PP
         DO 40 J = 1,NN
   40       BF(I+NN,J+MM) = C(I,J)
C
C     REDUCE THIS SYSTEM TO ONE WITH THE SAME INVARIANT ZEROS AND WITH
C     D OF FULL ROW RANK MU (THE NORMAL RANK OF THE ORIGINAL SYSTEM).
C
   45 RO = PP
      SIGMA = 0
      CALL RDUCE (BF,MAX,MM,NN,PP,EPS,RO,SIGMA,MU,NU,SUM,DUMMY)
      RANK = MU
      IF (NU .EQ. 0) RETURN
C
C     PERTRANSPOSE THE SYSTEM
C
      NUMU = NU+MU
      MNU = MM+NU
      NUMU1 = NUMU+1
      MNU1 = MNU+1
      DO 50 I = 1,NUMU
         DO 50 J = 1,MNU
   50       AF(MNU1-J,NUMU1-I) = BF(I,J)
      IF (MU .EQ. MM) GO TO 55
      PP = MM
      NN = NU
      MM = MU
C
C     REDUCE THE SYSTEM TO ONE WITH THE SAME INVARIANT ZEROS AND WITH
C     D SQUARE AND INVERTIBLE.
C
      RO = PP-MM
      SIGMA = MM
      CALL RDUCE (AF,MAX,MM,NN,PP,EPS,RO,SIGMA,MU,NU,SUM,DUMMY)
      IF (NU .EQ. 0) RETURN
C
C     PERFORM A UNITARY TRANSFORMATION ON THE COLUMNS OF
C        (LAMBDA*I - A  B)
C        (    -C        D)
C     IN ORDER TO REDUCE IT TO
C        (LAMBDA*BF - AF  X)
C        (     0          Y)
C     WITH Y AND BF SQUARE INVERTIBLE.
C
      MNU = MM+NU
   55 DO 70 I = 1,NU
         DO 60 J = 1,MNU
   60       BF(I,J) = 0.0D0
   70    BF(I,I+MM) = 1.0D0
      IF (RANK .EQ. 0) RETURN
      NU1 = NU+1
      I1 = NU+MU
      J1 = MNU+1
      I0 = MM
      DO 90 I = 1,MM
         I0 = I0-1
         DO 80 J = 1,NU1
   80       DUMMY(J) = AF(I1,I0+J)
      CALL HOUSH (DUMMY,NU1,NU1,EPS,ZERO,S)
      CALL TR2 (AF,MAX,DUMMY,S,1,I1,I0,NU1)
      CALL TR2 (BF,MAX,DUMMY,S,1,NU,I0,NU1)
   90 I1 = I1-1
      RETURN
      END