C
C      ________________________________________________________
C     |                                                        |
C     |      REDUCE A GENERAL MATRIX A TO BIDIAGONAL FORM      |
C     |     A = Q TIMES BIDIAGONAL MATRIX TIMES P TRANSPOSE    |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         D     --ARRAY WITH AT LEAST N ELEMENTS         |
C     |                                                        |
C     |         B     --ARRAY WITH AT LEAST M ELEMENTS         |
C     |                                                        |
C     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     |
C     |                                                        |
C     |         IQ    --AN INTEGER WHICH INDICATES WHICH COL-  |
C     |                 UMNS OF Q TO COMPUTE (= 0 MEANS NONE,  |
C     |                 = 1 MEANS FIRST L, = 2 MEANS LAST M-L, |
C     |                 = 3 MEANS ALL M WHERE L = MIN(M,N))    |
C     |                                                        |
C     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     |
C     |                                                        |
C     |         IP    --AN INTEGER (LIKE IQ) WHICH INDICATES   |
C     |                 WHICH COLUMNS OF P TO COMPUTE          |
C     |                                                        |
C     |         A     --ARRAY CONTAINING COEFFICIENT MATRIX    |
C     |                                                        |
C     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     |
C     |                                                        |
C     |         M     --ROW DIMENSION OF MATRIX STORED IN A    |
C     |                                                        |
C     |         N     --COLUMN DIMENSION OF MATRIX STORED IN A |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         D     --DIAGONAL OF BIDIAGONAL FORM            |
C     |                                                        |
C     |         B     --SUPERDIAGONAL (IF M .GE. N) OR SUBDIAG-|
C     |                 ONAL (IF M .LT. N) OF BIDIAGONAL FORM  |
C     |                                                        |
C     |         A     --THE HOUSEHOLDER VECTORS USED IN THE    |
C     |                 REDUCTION PROCESS                      |
C     |                                                        |
C     |         Q     --THE Q FACTOR IN THE BIDIAGONALIZATION  |
C     |                                                        |
C     |         P     --THE P FACTOR IN THE BIDIAGONALIZATION  |
C     |                                                        |
C     |    NOTE:                                               |
C     |                                                        |
C     |         EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT |
C     |         BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,|
C     |         THEN THE HOUSEHOLDER VECTORS IN A ARE DESTROYED|
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,MIN0,SQRT                    |
C     |    PACKAGE SUBROUTINES: HSR1-HSR5                      |
C     |________________________________________________________|
C
      SUBROUTINE BIDAG2(D,B,Q,LQ,IQ,P,LP,IP,A,LA,M,N)
      INTEGER H,I,IP,IQ,J,JP,JQ,K,L,LA,LP,LQ,M,N
      REAL A(LA,1),B(1),D(1),Q(LQ,1),P(LP,1),R,S,T,U
      L = MIN0(M,N)
      IF ( IQ .GE. 0 ) GOTO 20
10    WRITE(6,*) 'ERROR: INPUT PARAMETER IQ FOR SUBROUTINE BIDAG2'
      WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3'
      STOP
20    IF ( IQ .GT. 3 ) GOTO 10
      JQ = IQ
      IF ( IQ .LE. 1 ) GOTO 30
      IF ( IQ .EQ. 3 ) GOTO 30
      IF ( M .EQ. L ) JQ = 0
30    IF ( IP .GE. 0 ) GOTO 50
40    WRITE(6,*) 'ERROR: INPUT PARAMETER IP FOR SUBROUTINE BIDAG2'
      WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3'
      STOP
50    IF ( IP .GT. 3 ) GOTO 40
      JP = IP
      IF ( IP .LE. 1 ) GOTO 60
      IF ( IP .EQ. 3 ) GOTO 60
      IF ( N .EQ. L ) JP = 0
60    K = 1
      H = 2
      IF ( M .LT. N ) GOTO 330
      IF ( M .GT. 1 ) GOTO 70
      D(1) = A(1,1)
      IF ( IQ .GT. 0 ) Q(1,1) = 1.
      IF ( IP .GT. 0 ) P(1,1) = 1.
      RETURN
70    J = K
      K = H
      DO 80 I = K,M
80         IF ( A(I,J) .NE. 0. ) GOTO 110
      D(J) = A(J,J)
      A(J,J) = 0.
      DO 90 I = K,N
90         D(I) = A(J,I)
      IF ( JQ .EQ. 0 ) GOTO 200
      DO 100 I = J,M
100        Q(I,J) = 0.
      GOTO 200
110   T = ABS(A(J,J))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 130 L = I,M
           S = ABS(A(L,J))
           IF ( S .LE. T ) GOTO 120
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 130
120        R = R + (S*U)**2
130   CONTINUE
      S = T*SQRT(R)
      R = A(J,J)
      U = 1./SQRT(S*(S+ABS(R)))
      IF ( R .LT. 0. ) S = -S
      D(J) = -S
      A(J,J) = U*(R+S)
      DO 140 I = K,M
140        A(I,J) = A(I,J)*U
      IF ( JQ .EQ. 0 ) GOTO 160
      DO 150 I = J,M
150        Q(I,J) = A(I,J)
160   IF ( K .GT. N ) GOTO 620
      DO 190 L = K,N
           T = 0.
           DO 170 I = J,M
170             T = T + A(I,J)*A(I,L)
           A(J,L) = A(J,L) - T*A(J,J)
           D(L) = A(J,L)
           DO 180 I = K,M
180             A(I,L) = A(I,L) - T*A(I,J)
190   CONTINUE
200   H = K + 1
      IF ( K .LT. N ) GOTO 210
      IF ( K .GT. N ) GOTO 620
      IF ( M .EQ. N ) GOTO 610
      B(J) = A(J,N)
      GOTO 70
210   DO 220 I = H,N
220        IF ( D(I) .NE. 0. ) GOTO 240
      B(J) = D(K)
      A(J,K) = 0.
      IF ( IP .EQ. 0 ) GOTO 70
      DO 230 I = K,N
230        P(I,J) = 0.
      GOTO 70
240   T = ABS(D(K))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 260 L = I,N
           S = ABS(D(L))
           IF ( S .LE. T ) GOTO 250
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 260
250        R = R + (S*U)**2
260   CONTINUE
      S = T*SQRT(R)
      R = D(K)
      U = 1./SQRT(S*(S+ABS(R)))
      IF ( R .LT. 0. ) S = -S
      D(K) = U*(R+S)
      DO 270 I = H,N
270        D(I) = D(I)*U
      IF ( IP .EQ. 0 ) GOTO 290
      DO 280 I = K,N
280        P(I,J) = D(I)
290   B(J) = -S
      DO 300 I = K,M
300        B(I) = 0.
      DO 310 L = K,N
           T = D(L)
           A(J,L) = T
           DO 310 I = K,M
310             B(I) = B(I) + T*A(I,L)
      DO 320 L = K,N
           T = D(L)
           DO 320 I = K,M
320             A(I,L) = A(I,L) - T*B(I)
      GOTO 70
330   DO 340 I = K,N
340        D(I) = A(K,I)
350   J = K
      K = H
      DO 360 I = K,N
360        IF ( D(I) .NE. 0. ) GOTO 370
      U = D(J)
      D(J) = 0.
      GOTO 440
370   T = ABS(D(J))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 390 L = I,N
           S = ABS(D(L))
           IF ( S .LE. T ) GOTO 380
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 390
380        R = R + (S*U)**2
390   CONTINUE
      S = T*SQRT(R)
      R = D(J)
      U = 1./SQRT(S*(S+ABS(R)))
      IF ( R .LT. 0. ) S = -S
      D(J) = U*(R+S)
      DO 400 I = K,N
400        D(I) = D(I)*U
      U = -S
      IF ( K .GT. M ) GOTO 470
      DO 410 I = K,M
410        B(I) = 0.
      DO 420 L = J,N
           T = D(L)
           A(J,L) = T
           DO 420 I = K,M
420             B(I) = B(I) + T*A(I,L)
      DO 430 L = J,N
           T = D(L)
           DO 430 I = K,M
430             A(I,L) = A(I,L) - T*B(I)
440   H = K + 1
      IF ( IP .EQ. 0 ) GOTO 460
      DO 450 I = J,N
450        P(I,J) = D(I)
460   D(J) = U
      IF ( K .LT. M ) GOTO 490
      IF ( K .GT. M ) GOTO 620
      B(J) = A(M,J)
      GOTO 330
470   DO 480 I = J,N
480        A(J,I) = D(I)
      GOTO 440
490   DO 500 I = H,M
500        IF ( A(I,J) .NE. 0. ) GOTO 520
      B(J) = A(K,J)
      A(K,J) = 0.
      IF ( IQ .EQ. 0 ) GOTO 330
      DO 510 I = K,M
510        Q(I,J) = 0.
      GOTO 330
520   T = ABS(A(K,J))
      IF ( T .NE. 0. ) U = 1./T
      R = 1.
      DO 540 L = I,M
           S = ABS(A(L,J))
           IF ( S .LE. T ) GOTO 530
           U = 1./S
           R = 1. + R*(T*U)**2
           T = S
           GOTO 540
530        R = R + (S*U)**2
540   CONTINUE
      S = T*SQRT(R)
      R = A(K,J)
      U = 1./SQRT(S*(S+ABS(R)))
      IF ( R .LT. 0. ) S = -S
      B(J) = -S
      A(K,J) = U*(R+S)
      DO 550 I = H,M
550        A(I,J) = A(I,J)*U
      IF ( IQ .EQ. 0 ) GOTO 570
      DO 560 I = K,M
560        Q(I,J) = A(I,J)
570   DO 600 L = K,N
           T = 0.
           DO 580 I = K,M
580             T = T + A(I,J)*A(I,L)
           A(K,L) = A(K,L) - T*A(K,J)
           D(L) = A(K,L)
           DO 590 I = H,M
590             A(I,L) = A(I,L) - T*A(I,J)
600   CONTINUE
      GOTO 350
610   D(N) = A(N,N)
      B(N-1) = A(N-1,N)
620   IF ( JQ .EQ. 0 ) GOTO 650
      IF ( N .GT. M ) GOTO 640
      IF ( N .EQ. M ) GOTO 630
      IF ( JQ .EQ. 1 ) CALL HSR3(Q,LQ,M,N)
      IF ( JQ .EQ. 2 ) CALL HSR4(Q,LQ,M,N)
      IF ( JQ .EQ. 3 ) CALL HSR5(Q,LQ,M,N)
      GOTO 650
630   CALL HSR2(Q,LQ,M)
      GOTO 650
640   CALL HSR1(Q,LQ,M)
650   IF ( JP .EQ. 0 ) RETURN
      IF ( N .LE. M ) GOTO 660
      IF ( JP .EQ. 1 ) CALL HSR3(P,LP,N,M)
      IF ( JP .EQ. 2 ) CALL HSR4(P,LP,N,M)
      IF ( JP .EQ. 3 ) CALL HSR5(P,LP,N,M)
      RETURN
660   CALL HSR1(P,LP,N)
      RETURN
      END