SUBROUTINE RDUCE (ABCD,MDIMA,M,N,P,EPS,RO,SIGMA,MU,NU,SUM,
     +                   DUMMY)
C
C        THIS ROUTINE EXTRACTS FROM THE (N+P)X(N+M) SYSTEM (B A)
C                                                          (D C)
C        A (NU+MU)X(M+NU) 'REDUCED' SYSTEM (B' A') HAVING THE SAME
C                                          (D' C')
C        TRANSMISSION ZEROS BUT WITH D' OF FULL ROW RANK.  THE
C        SYSTEM (A',B',C',D') OVERWRITES THE OLD SYSTEM.  EPS IS
C        THE NOISE LEVEL.  SUM(MAX(M,P)) AND DUMMY(MAX(N,P)) (SIC)
C        ARE WORKING ARRAYS.
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)
      INTEGER TAU,P,RO,RO1,SIGMA
      LOGICAL ZERO
      DIMENSION ABCD(MDIMA,1),DUMMY(1),SUM(1)
      MU = P
      NU = N
   10 IF (MU .EQ. 0) RETURN
      RO1 = RO
      MNU = M+NU
      NUMU = NU+MU
      IF (M .EQ. 0) GO TO 120
      RO1 = RO1+1
      IROW = NU
      IF (SIGMA .LE. 1) GO TO 40
C
C     COMPRESS ROWS OF D.  FIRST EXPLOIT TRIANGULAR SHAPE.
C
      M1 = SIGMA-1
      DO 30 ICOL = 1,M1
         DO 20 J = 1,RO1
   20       DUMMY(J) = ABCD(IROW+J,ICOL)
         CALL HOUSH (DUMMY,RO1,1,EPS,ZERO,S)
         IICOL = ICOL
         CALL TR1 (ABCD,MDIMA,DUMMY,S,IROW,RO1,IICOL,MNU)
   30    IROW = IROW+1
C
C     CONTINUE WITH HOUSEHOLDER WITH PIVOTING
C
   40 IF (SIGMA .NE. 0) GO TO 45
      SIGMA = 1
      RO1 = RO1-1
   45 IF (SIGMA .EQ. M) GO TO 60
      DO 55 ICOL = SIGMA,M
         DUM = 0.0D0
         DO 50 J = 1,RO1
   50       DUM = DUM+ABCD(IROW+J,ICOL)*ABCD(IROW+J,ICOL)
   55    SUM(ICOL) = DUM
   60 DO 100 ICOL = SIGMA,M
C
C        PIVOT IF NECESSARY
C
         IICOL = ICOL
         IF (ICOL .EQ. M) GO TO 80
         CALL PIVOT (SUM,DUM,IBAR,IICOL,M)
         IF (IBAR .EQ. ICOL) GO TO 80
         SUM(IBAR) = SUM(ICOL)
         SUM(ICOL) = DUM
         DO 70 I = 1,NUMU
            DUM = ABCD(I,ICOL)
            ABCD(I,ICOL) = ABCD(I,IBAR)
   70       ABCD(I,IBAR) = DUM
C
C        PERFORM HOUSEHOLDER TRANSFORMATION
C
   80    DO 90 I = 1,RO1
   90       DUMMY(I) = ABCD(IROW+I,ICOL)
         CALL HOUSH (DUMMY,RO1,1,EPS,ZERO,S)
         IF (ZERO) GO TO 120
         IF (RO1 .EQ. 1) RETURN
         CALL TR1 (ABCD,MDIMA,DUMMY,S,IROW,RO1,IICOL,MNU)
         IROW = IROW+1
         RO1 = RO1-1
         DO 100 J = ICOL,M
  100       SUM(J) = SUM(J)-ABCD(IROW,J)*ABCD(IROW,J)
  120 TAU = RO1
      SIGMA = MU-TAU
C
C     COMPRESS THE COLUMNS OF C
C
      IF (NU .LE. 0) GO TO 220
      I1 = NU+SIGMA
      MM1 = M+1
      N1 = NU
      IF (TAU .EQ. 1) GO TO 140
      DO 135 I = 1,TAU
         DUM = 0.0D0
         DO 130 J = MM1,MNU
  130       DUM = DUM+ABCD(I1+I,J)*ABCD(I1+I,J)
  135    SUM(I) = DUM
  140 DO 200 RO1 = 1,TAU
         RO = RO1-1
         I = TAU-RO
         I2 = I+I1
C
C        PIVOT IF NECESSARY
C
         IF (I .EQ. 1) GO TO 160
         CALL PIVOT (SUM,DUM,IBAR,1,I)
         IF (IBAR .EQ. I) GO TO 160
         SUM(IBAR) = SUM(I)
         SUM(I) = DUM
         DO 150 J = MM1,MNU
            DUM = ABCD(I2,J)
            ABCD(I2,J) = ABCD(IBAR+I1,J)
  150       ABCD(IBAR+I1,J) = DUM
C
C        PERFORM HOUSEHOLDER TRANSFORMATION
C
  160    DO 170 J = 1,N1
  170       DUMMY(J) = ABCD(I2,M+J)
         CALL HOUSH (DUMMY,N1,N1,EPS,ZERO,S)
         IF (ZERO) GO TO 210
         IF (N1 .EQ. 1) GO TO 200
         CALL TR2 (ABCD,MDIMA,DUMMY,S,1,I2,M,N1)
         MN1 = M+N1
         CALL TR1 (ABCD,MDIMA,DUMMY,S,0,N1,1,MN1)
         DO 190 J = 1,I
  190       SUM(J) = SUM(J)-ABCD(I1+J,MN1)*ABCD(I1+J,MN1)
         MNU = MNU-1
  200    N1 = N1-1
      RO = TAU
  210 NU = NU-RO
      MU = SIGMA+RO
      IF (RO .EQ. 0) RETURN
      GO TO 10
  220 MU = SIGMA
      NU = 0
      RETURN
      END