SUBROUTINE ZM1(NA,A,NB,B,NC,C)
C  Multiply two z-transforms, C(z)=A(z)*B(z)
C  Inputs:
C  NA    Highest power of A
C  A     Coefficients A(0) through A(NA)
C  NB    Hightest power of B
C  B     Coefficients  B(0) through B(NB)
C  NC    Highest power of C to compute
C  Outputs:
C  C     Coefficients C(0) through C(NC)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(0:NA),B(0:NB),C(0:NC)
      MC=MIN0(NA+NB,NC)
      IF(MC.LT.NC)CALL DSET(NC-MC,0.0D0,C(MC+1),1)
      DO 20 IC=0,MC
      AA=0.0D0
      K1=MAX0(0,IC-NB)
      K2=MIN0(IC,NA)
      DO 10 K=K1,K2
   10 AA=AA+A(K)*B(IC-K)
   20 C(IC)=AA
      RETURN
      END

      SUBROUTINE ZQ1(NA,A,NB,B,NC,C)
C  Quotient of two z-transforms, C(z)=A(z)/B(z)
C  Inputs:
C  NA    Highest power of A
C  A     Coefficients A(0) through A(NA)
C  NB    Hightest power of B
C  B     Coefficients  B(0) through B(NB)
C  NC    Highest power of C to compute
C  Outputs:
C  C     Coefficients C(0) through C(NC)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(0:NA),B(0:NB),C(0:NC)
      C(0)=A(0)/B(0)
      DO 20 IC=1,NC
      AA=0.0D0
      IF(IC.LE.NA)AA=A(IC)
      K2=MIN0(IC,NB)
      DO 10 K=1,K2
   10 AA=AA-B(K)*C(IC-K)
   20 C(IC)=AA/B(0)
      RETURN
      END

      SUBROUTINE ZI1(NA,A,NB,B)
C  Inverse of z-transform, B(z)=1/A(z)
C  Inputs:
C  NA    Highest power of A
C  A     Coefficients A(0) through A(NA)
C  NB    Highest power of B to compute
C  Outputs:
C  B     Coefficients B(0) through B(NB)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(0:NA),B(0:NB)
      B(0)=1.0D0/A(0)
      DO 20 IB=1,NB
      AA=0.0D0
      K2=MIN0(IB,NA)
      DO 10 K=1,K2
   10 AA=AA+A(K)*B(IB-K)
   20 B(IB)=-B(0)*AA
      RETURN
      END

      SUBROUTINE ZM2(NA,NRA,NCA,A,LDRA,LDCA,NB,NRB,NCB,B,LDRB,LDCB,
     1               NC,NRC,NCC,C,LDRC,LDCC)
C  Multiply two matrix z-transforms, C(z)=A(z)*B(z)
C  Inputs:
C  NA    Highest power of A
C  NRA   Number of rows in A
C  NCA   Number of columns in A
C  A     Coefficients A(0) through A(NA)
C  LDRA  Leading row dimension of A
C  LDCA  Leading column dimension of A
C  NB    Hightest power of B
C  NRB   Number of rows in B
C  NCB   Number of columns in B
C  B     Coefficients  B(0) through B(NB)
C  LDRB  Leading row dimension of B
C  LDCB  Leading column dimension of B
C  NC    Highest power of C to compute
C  NRC   Number of rows in C
C  NCC   Number of columns in C
C  LDRC  Leading row dimension of C
C  LDCC  leading column dimension of C
C  Outputs:
C  C     Coefficients C(0) through C(NC)
      PARAMETER(LD=100)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB),C(LDRC,LDCC,0:NC)
      MC=MIN0(NA+NB,NC)
      DO 10 IC=0,NC
   10 CALL UM0SET(NRC,NCC,C(1,1,IC),LDRC)
      DO 20 IC=0,MC
      K1=MAX0(0,IC-NB)
      K2=MIN0(IC,NA)
      DO 20 K=K1,K2
      CALL DMRRRR(NRA,NRB,A(1,1,K),LDRA,NRB,NCB,B(1,1,IC-K),LDRB,
     1            NRC,NCC,A1,LD)
   20 CALL UMADD(NRC,NCC,C(1,1,IC),LDRC,A1,LD,C(1,1,IC),LDRC)
      RETURN
      END

      SUBROUTINE ZQ2(NA,NRA,NCA,A,LDRA,LDCA,NB,NRB,NCB,B,LDRB,LDCB,
     1               NC,NRC,NCC,C,LDRC,LDCC)
C  Quotient of two matrix z-transforms, C(z)=[B(z)^-1]*A(z)
C  Inputs:
C  NA    Highest power of A
C  NRA   Number of rows in A
C  NCA   Number of columns in A
C  A     Coefficients A(0) through A(NA)
C  LDRA  Leading row dimension of A
C  LDCA  Leading column dimension of A
C  NB    Hightest power of B
C  NRB   Number of rows in B
C  NCB   Number of columns in B
C  B     Coefficients  B(0) through B(NB)
C  LDRB  Leading row dimension of B
C  LDCB  Leading column dimension of B
C  NC    Highest power of C to compute
C  NRC   Number of rows in C
C  NCC   Number of columns in C
C  LDRC  Leading row dimension of C
C  LDCC  leading column dimension of C
C  Outputs:
C  C     Coefficients C(0) through C(NC)
      PARAMETER(LD=100)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB),C(LDRC,LDCC,0:NC)
      CALL DLINRG(NRB,B,LDRB,A2,LD)
      CALL DMRRRR(NRB,NCB,A2,LD,NRA,NCA,A,LDRA,NRC,NCC,C,LDRC)
      DO 20 IC=1,NC
      CALL UM0SET(NRA,NCA,A1,LD)
      IF(IC.LE.NA)CALL UCRGRG(NRA,NCA,A(1,1,IC),LDRA,A1,LD)
      K2=MIN0(IC,NB)
      DO 10 K=1,K2
      CALL DMRRRR(NRB,NCB,B(1,1,K),LDRB,NRC,NCC,C(1,1,IC-K),LDRC,
     1            NRC,NCC,C(1,1,IC),LDRC)
   10 CALL UMSUB(NRC,NCC,A1,LD,C(1,1,IC),LDRC,A1,LD)
   20 CALL DMRRRR(NRB,NCB,A2,LD,NRC,NCC,A1,LD,NRC,NCC,C(1,1,IC),LDRC)
      RETURN
      END

      SUBROUTINE ZI2(M,NA,A,LDRA,LDCA,NB,B,LDRB,LDCB)
C  Invert a matrix z-transform, B(z)=A(z)^-1
C  Inputs:
C  NA    Highest power of A
C  NRA   Number of rows in A
C  NCA   Number of columns in A
C  A     Coefficients A(0) through A(NA)
C  LDRA  Leading row dimension of A
C  LDCA  Leading column dimension of A
C  NB    Hightest power of B to compute
C  NRB   Number of rows in B
C  NCB   Number of columns in B
C  LDRB  Leading row dimension of B
C  LDCB  Leading column dimension of B
C  Outputs:
C  B     Coefficients B(0) through B(NB)
      PARAMETER(LD=100)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION A(LDRA,LDCA,0:NA),B(LDRB,LDCB,0:NB)
      CALL DLINRG(M,A,LDRA,B,LDRB)
      DO 20 IB=1,NB
      CALL UM0SET(M,M,A1,LD)
      K2=MIN0(IB,NA)
      DO 10 K=1,K2
      CALL DMRRRR(M,M,A(1,1,K),LDRA,M,M,B(1,1,IB-K),LDRB,
     1            M,M,B(1,1,IB),LDRB)
   10 CALL UMSUB(M,M,A1,LD,B(1,1,IB),LDRB,A1,LD)
   20 CALL DMRRRR(M,M,B,LDRB,M,M,A1,LD,M,M,B(1,1,IB),LDRB)
      RETURN
      END