C
C      ________________________________________________________
C     |                                                        |
C     |  COMPUTE SINGULAR VALUE DECOMPOSITION OF A BIDIAGONAL  |
C     | MATRIX:  A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         D     --DIAGONAL                               |
C     |                                                        |
C     |         N     --DIMENSION OF BIDIAGONAL MATRIX         |
C     |                                                        |
C     |         U     --OFF DIAGONAL                           |
C     |                                                        |
C     |         IU    --= 0 IF U IS SUPERDIAGONAL AND = 1 IF   |
C     |                 U IS SUBDIAGONAL                       |
C     |                                                        |
C     |         Q     --EITHER A SEGMENT OF AN IDENTITY MATRIX |
C     |                 OR A SEGMENT OF THE MATRIX USED TO     |
C     |                 BIDIAGONALIZE THE COEFFICIENT MATRIX   |
C     |                                                        |
C     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     |
C     |                                                        |
C     |         MQ    --NUMBER OF MATRIX ELEMENTS CONTAINED IN |
C     |                 EACH COLUMN OF Q                       |
C     |                                                        |
C     |         IQ    --AN INTEGER WHICH INDICATES WHETHER     |
C     |                 COLUMNS OF Q TO BE PROCESSED (= 0 MEANS|
C     |                 NO AND = 1 MEANS YES)                  |
C     |                                                        |
C     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     |
C     |                                                        |
C     |         MP    --NUMBER OF MATRIX ELEMENTS CONTAINED IN |
C     |                 EACH COLUMN OF P                       |
C     |                                                        |
C     |         IP    --AN INTEGER DEFINED LIKE IQ             |
C     |                                                        |
C     |         E     --WORK ARRAY WITH AT LEAST N ELEMENTS    |
C     |                                                        |
C     |         F     --WORK ARRAY WITH AT LEAST N ELEMENTS    |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         Q     --Q FACTOR IN THE SINGULAR VALUE DECOMP. |
C     |                                                        |
C     |         D     --SINGULAR VALUES IN DECREASING ORDER    |
C     |                                                        |
C     |         P     --P FACTOR IN THE SINGULAR VALUE DECOMP. |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,AMAX1,SIGN,SQRT              |
C     |    PACKAGE SUBROUTINES: EIG3,FGV,HSR1-HSR5,SCL,SFT,    |
C     |                         SNG0,SNG1,SORT2                |
C     |________________________________________________________|
C
      SUBROUTINE SINGB(D,N,U,IU,Q,LQ,MQ,IQ,P,LP,MP,IP,E,F)
      INTEGER G,H,I,ID,IP,IQ,IU,J,JL,JR,J0,K,K2,L,LL,LP,LQ,L1
      INTEGER M,MP,MQ,N,NS
      REAL D(1),E(1),F(1),Q(LQ,1),P(LP,1),U(1)
      REAL B,C,R,S,T,T0,T1,T2,T3,V,W,X,Y,Z
      IF ( N .GT. 1 ) GOTO 10
      IF ( IQ .EQ. 1 ) Q(1,1) = 1.
      IF ( IP .EQ. 1 ) P(1,1) = 1.
      IF ( D(1) .GE. 0. ) RETURN
      D(1) = -D(1)
      IF ( IQ .EQ. 1 ) Q(1,1) = -1.
      RETURN
10    JL = IQ
      IF ( JL .EQ. 0 ) JL = 3
      IF ( JL .NE. 3 ) JL = 1
      JR = IP
      IF ( JR .EQ. 0 ) JR = 3
      IF ( JR .NE. 3 ) JR = 2
      IF ( IU .EQ. 0 ) GOTO 20
      I = JL
      JL = JR
      JR = I
20    J = 0
      L = N - 1
      K2 = N - 2
      DO 30 I = 1,L
           E(I) = 1.
           F(I) = 1.
           IF ( U(I) .EQ. 0. ) J = I
30         IF ( D(I) .EQ. 0. ) J = I
      E(N) = 1.
      F(N) = 1.
      B = 65536.**(-3)
      T = 1.
40    T = .5*T
      S = 1. + T
      IF ( S .GT. 1. ) GOTO 40
      T0 = 1./(T+T)
      T2 = (T+T)**2
      NS = 50*N
      L1 = 0
      K = N
      LL = 0
      GOTO 70
50    J = 0
      DO 60 I = 1,L
           IF ( U(I) .EQ. 0. ) J = I
60         IF ( D(I) .EQ. 0. ) J = I
70    IF ( J .EQ. 0 ) GOTO 140
      IF ( U(J) .EQ. 0. ) GOTO 140
C     -------------------------------
C     |*** ZERO DIAGONAL ELEMENT ***|
C     -------------------------------
      I = J
      V = U(J)
      U(J) = 0.
      S = -V*ABS(E(J))
C     ---------------------------
C     |*** PROCESS LEFT SIDE ***|
C     ---------------------------
80    H = I
      I = I + 1
      R = ABS(E(I))*D(I)
      CALL FGV(X,Y,T,R,S,E(J),E(I))
      IF ( T .EQ. 1. ) GOTO 90
      D(I) = D(I) - Y*V
      CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,J,I,X,Y)
      IF ( I .EQ. K ) GOTO 100
      V = X*U(I)
      S = -V*ABS(E(J))
      GOTO 80
90    D(I) = D(I)*Y - V
      CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,J,I,X,Y)
      IF ( I .EQ. K ) GOTO 100
      V = U(I)
      U(I) = V*Y
      S = -V*ABS(E(J))
      GOTO 80
100   IF ( J .EQ. 1 ) GOTO 130
      I = J - 1
      S = U(I)
      U(I) = 0.
C     ----------------------------
C     |*** PROCESS RIGHT SIDE ***|
C     ----------------------------
110   H = I
      I = I - 1
      R = D(H)
      CALL FGV(X,Y,T,R,S,F(H),F(J))
      IF ( T .EQ. 1. ) GOTO 120
      D(H) = R + X*S
      CALL SNG0(Q,LQ,MQ,P,LP,MP,JR,H,J,X,Y)
      IF ( H .EQ. 1 ) GOTO 130
      S = -Y*U(I)
      GOTO 110
120   D(H) = X*R + S
      CALL SNG1(Q,LQ,MQ,P,LP,MP,JR,H,J,X,Y)
      IF ( H .EQ. 1 ) GOTO 130
      S = -U(I)
      U(I) = X*U(I)
      GOTO 110
130   CALL SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,1,K,JL,JR)
140   J = J + 1
      IF ( J .EQ. K ) GOTO 320
      S = 0.
      T = 0.
C     -----------------------------
C     |*** SET ERROR TOLERANCE ***|
C     -----------------------------
      DO 150 I = J,L
           X = ABS((D(I)*E(I))*(D(I)*F(I)))
           Y = ABS((U(I)*E(I))*(U(I)*F(I+1)))
           S = AMAX1(S,X,Y)
150        T = T + X + Y
      X = ABS((E(K)*D(K))*(F(K)*D(K)))
      S = AMAX1(S,X)
      T3 = S*T2
      T = T + X
      IF ( T .EQ. 0. ) T = 1.
      T1 = T0/T
      GOTO 280
160   LL = LL + 1
      IF ( LL .GT. NS ) GOTO 530
      IF ( L .GT. J ) GOTO 170
      S = 0.
      T = 0.
      GOTO 180
170   S = U(K2)
      T = E(K2)
C     -----------------------
C     |*** COMPUTE SHIFT ***|
C     -----------------------
180   CALL SFT(C,D(K),D(L),U(L),S,E(K),E(L),T,F(K),F(L))
      V = ABS(E(J))
      W = ABS(F(J))
      T = D(J)*W
      R = T*(D(J)*V) - C
      S = T*(U(J)*V)
      ID = 1
      J0 = J
      I = J
      H = J - 1
190   G = H
      H = I
      I = I + 1
      Z = ABS(F(I))
C     ----------------------------
C     |*** PROCESS RIGHT SIDE ***|
C     ----------------------------
      CALL FGV(X,Y,T,R,S,F(H),F(I))
      V = D(H)
      IF ( T .EQ. 1 ) GOTO 210
      IF ( H .EQ. J ) GOTO 200
      T = U(G) + X*S
      U(G) = T
      IF ( ABS((T*E(G))*(T*F(H))) .GT. T3 ) GOTO 200
      J = H
      ID = 1
200   R = V + X*U(H)
      S = X*D(I)
      U(H) = U(H) - Y*V
      CALL SNG0(Q,LQ,MQ,P,LP,MP,JR,H,I,X,Y)
      GOTO 230
210   IF ( H .EQ. J ) GOTO 220
      T = X*U(G) + S
      U(G) = T
      IF ( ABS((T*E(G))*(T*F(H))) .GT. T3 ) GOTO 220
      J = H
      ID = 1
220   R = X*V + U(H)
      S = D(I)
      U(H) = Y*U(H) - V
      D(I) = Y*S
      CALL SNG1(Q,LQ,MQ,P,LP,MP,JR,H,I,X,Y)
C     ---------------------------
C     |*** PROCESS LEFT SIDE ***|
C     ---------------------------
230   CALL FGV(X,Y,T,R,S,E(H),E(I))
      IF ( T .EQ. 1. ) GOTO 250
      T = R + X*S
      D(H) = T
      IF ( ABS((T*E(H))*(T*F(H))) .GT. T3 ) GOTO 240
      ID = 0
      J = H
240   R = U(H) + X*D(I)
      D(I) = D(I) - Y*U(H)
      U(H) = R
      CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,H,I,X,Y)
      IF ( I .EQ. K ) GOTO 270
      S = X*U(I)
      GOTO 190
250   T = S + X*R
      D(H) = T
      IF ( ABS((T*E(H))*(T*F(H))) .GT. T3 ) GOTO 260
      ID = 0
      J = H
260   R = D(I) + X*U(H)
      D(I) = Y*D(I) - U(H)
      U(H) = R
      CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,H,I,X,Y)
      IF ( I .EQ. K ) GOTO 270
      S = U(I)
      U(I) = S*Y
      GOTO 190
270   CALL SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,J0,K,JL,JR)
      IF ( ID .EQ. 0 ) GOTO 70
280   W = E(L)
      X = E(K)
      Y = F(L)
      Z = F(K)
      R = ABS((X*D(K))*(Z*D(K)))
      S = ABS((W*D(L))*(Y*D(L)))
      T = ABS((X*U(L))*(Y*U(L)))
C     ------------------------------
C     |*** TEST FOR CONVERGENCE ***|
C     ------------------------------
      IF ( (S*T1)*(T*T1) .GT. 1. ) GOTO 160
      L1 = L1 + 1
      IF ( L1 .GT. 40 ) GOTO 290
      IF ( T .EQ. 0. ) GOTO 290
      R = R + T
      IF ( (S/R)*(T/R) .GT. T2 ) GOTO 160
290   L1 = 0
      IF ( S .GT. R ) GOTO 310
      R = -D(K)*ABS(X)
      S = U(L)*ABS(W)
      CALL FGV(W,Y,T,R,S,E(L),E(K))
      X = E(K)
      IF ( T .EQ. 1. ) GOTO 300
      D(K) = D(K) - Y*U(L)
      CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,L,K,W,Y)
      GOTO 310
300   D(L) = D(L)*W
      D(K) = Y*D(K) - U(L)
      CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,L,K,W,Y)
310   T = SIGN(SQRT(ABS(X)),X)*D(K)*SIGN(SQRT(ABS(Z)),Z)
      IF ( T .LT. 0. ) E(K) = -E(K)
      D(K) = ABS(T)
      K = L
      L = K2
      K2 = K2 - 1
      IF ( K .GT. J ) GOTO 280
320   X = E(K)
      Z = F(K)
      T = SIGN(SQRT(ABS(X)),X)*D(K)*SIGN(SQRT(ABS(Z)),Z)
      IF ( T .LT. 0. ) E(K) = -E(K)
      D(K) = ABS(T)
      K = L
      L = K2
      K2 = K2 - 1
      IF ( K .GT. 1 ) GOTO 50
      IF ( K .EQ. 0 ) GOTO 330
      GOTO 320
C     -------------------------
C     |*** RESCALE P AND Q ***|
C     -------------------------
330   GOTO (340,360,380),JL
340   DO 350 J = 1,N
           T = E(J)
           T = SIGN(SQRT(ABS(T)),T)
           DO 350 I = 1,MQ
350             Q(I,J) = Q(I,J)*T
      GOTO 380
360   DO 370 J = 1,N
           T = E(J)
           T = SIGN(SQRT(ABS(T)),T)
           DO 370 I = 1,MP
370             P(I,J) = P(I,J)*T
380   GOTO (390,410,430),JR
390   DO 400 J = 1,N
           T = F(J)
           T = SIGN(SQRT(ABS(T)),T)
           DO 400 I = 1,MQ
400             Q(I,J) = Q(I,J)*T
      GOTO 430
410   DO 420 J = 1,N
           T = F(J)
           T = SIGN(SQRT(ABS(T)),T)
           DO 420 I = 1,MP
420             P(I,J) = P(I,J)*T
C     --------------------------------
C     |*** REORDER THE EIGENPAIRS ***|
C     --------------------------------
430   CALL SORT2(D,E,F,N)
      DO 440 I = 1,N
           J = E(I)
440        F(J) = I
      M = N + 1
      DO 520 J = 1,N
           L = M - J
           K = E(L)
           I = F(J)
           E(I) = K
           F(K) = I
           T = D(J)
           D(J) = D(K)
           D(K) = T
           IF ( JL .EQ. 1 ) GOTO 450
           IF ( JR .NE. 1 ) GOTO 480
450        S = 0.
           DO 460 I = 1,MQ
                T = Q(I,K)
                Q(I,K) = Q(I,J)
                Q(I,J) = T
460             S = S + T*T
           S = 1./SQRT(S)
           DO 470 I = 1,MQ
470             Q(I,J) = S*Q(I,J)
480        IF ( JR .EQ. 2 ) GOTO 490
           IF ( JL .NE. 2 ) GOTO 520
490        S = 0.
           DO 500 I = 1,MP
                T = P(I,K)
                P(I,K) = P(I,J)
                P(I,J) = T
500             S = S + T*T
           S = 1./SQRT(S)
           DO 510 I = 1,MP
510             P(I,J) = S*P(I,J)
520   CONTINUE
      E(1) = N
      RETURN
530   K = N - K + 1
      WRITE(6,*) 'SINCE THE STOPPING CRITERION NOT SATISFIED'
      WRITE(6,*) 'AFTER',NS,'ITERATIONS, WE STOP WHILE COMPUTING'
      WRITE(6,*) 'EIGENVALUE NUMBER',K
      E(1) = K
      RETURN
      END
C%
      SUBROUTINE FGV(X,Y,S,P,Q,A,B)
      REAL A,B,C,P,Q,R,S,T,X,Y
      IF ( ABS(P) .GT. ABS(Q) ) GOTO 10
      IF ( Q .EQ. 0. ) GOTO 110
      R = A/B
      S = P/Q
      T = ABS(R)*S*S
      IF ( T .LT. 1 ) GOTO 70
      T = T/(1.+T)
      R = B/A
      S = Q/P
      GOTO 20
10    R = B/A
      S = Q/P
      T = ABS(R)*S*S
      IF ( T .GT. 1. ) GOTO 60
      T = 1./(1.+T)
20    A = SIGN(A*T,P)
      IF ( SIGN(R,P) .EQ. R ) GOTO 40
      B = -ABS(B*T)
      GOTO 50
40    B = ABS(B*T)
50    Y = S
      X = ABS(R)*S
      S = 0.
      RETURN
60    T = T/(1.+T)
      R = A/B
      S = P/Q
      GOTO 80
70    T = 1./(1.+T)
80    C = A
      A = SIGN(B*T,Q)
      IF ( SIGN(R,Q) .EQ. R ) GOTO 90
      B = -ABS(C*T)
      GOTO 100
90    B = ABS(C*T)
100   Y = S
      X = ABS(R)*S
      S = 1.
      RETURN
110   X = 0.
      Y = 0.
      S = 0.
      RETURN
      END
C%
      SUBROUTINE SNG0(Q,LQ,M,P,LP,N,L,J,K,X,Y)
      INTEGER I,J,K,L,LP,LQ,M,N
      REAL Q(LQ,1),P(LP,1),S,T,X,Y
      GOTO (10,30,50),L
10    DO 20 I = 1,M
           T = Q(I,J)
           S = Q(I,K)
           Q(I,J) = T + X*S
20         Q(I,K) = S - Y*T
      RETURN
30    DO 40 I = 1,N
           T = P(I,J)
           S = P(I,K)
           P(I,J) = T + X*S
40         P(I,K) = S - Y*T
50    RETURN
      END
C%
      SUBROUTINE SNG1(Q,LQ,M,P,LP,N,L,J,K,X,Y)
      INTEGER I,J,K,L,LP,LQ,M,N
      REAL Q(LQ,1),P(LP,1),S,T,X,Y
      GOTO (10,30,50),L
10    DO 20 I = 1,M
           T = Q(I,J)
           S = Q(I,K)
           Q(I,J) = X*T + S
20         Q(I,K) = Y*S - T
      RETURN
30    DO 40 I = 1,N
           T = P(I,J)
           S = P(I,K)
           P(I,J) = X*T + S
40         P(I,K) = Y*S - T
50    RETURN
      END
C%
      SUBROUTINE SFT(S,A,B,C,D,E2,E1,E0,F2,F1)
      REAL A,B,C,D,E0,E1,E2,F1,F2,G0,G1,G2,H1,H2,S,W,X,Y,Z
      G0 = ABS(E0)
      G1 = ABS(E1)
      G2 = ABS(E2)
      H1 = ABS(F1)
      H2 = ABS(F2)
      W = (A*G2)*(A*H2) + (C*G1)*(C*H2)
      X = (B*G1)*(B*H1) + (D*G0)*(D*H1)
      Y = (B*G1)*(C*H1)
      Z = (B*G1)*(C*H2)
      CALL EIG3(S,S,X,W,Y,Z)
      RETURN
      END
C%
      SUBROUTINE SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,J,K,JL,JR)
      INTEGER I,J,JL,JR,K,L,LP,LQ,M,MP,MQ,N
      REAL D(1),E(1),F(1),P(LP,1),Q(LQ,1),U(1),B,R,S,T,V
      INTEGER H
      T = 1.
      DO 10 I = J,K
           IF ( ABS(E(I)) .LT. T ) T = ABS(E(I))
10         IF ( ABS(F(I)) .LT. T ) T = ABS(F(I))
      IF ( T .GT. B ) RETURN
C     ----------------------------
C     |*** RESCALE THE MATRIX ***|
C     ----------------------------
      R = E(J)
      R = SIGN(SQRT(ABS(R)),R)
      E(J) = 1.
      S = F(J)
      S = SIGN(SQRT(ABS(S)),S)
      F(J) = 1.
      D(J) = R*D(J)*S
      T = R
      IF ( JL .EQ. 1 ) GOTO 20
      IF ( JR .NE. 1 ) GOTO 40
      T = S
20    DO 30 I = 1,MQ
30         Q(I,J) = Q(I,J)*T
40    T = S
      IF ( JR .EQ. 2 ) GOTO 50
      IF ( JL .NE. 2 ) GOTO 70
      T = R
50    DO 60 I = 1,MP
60         P(I,J) = P(I,J)*T
70    L = J + 1
      H = J
      DO 130 M = L,K
           T = E(M)
           T = SIGN(SQRT(ABS(T)),T)
           E(M) = 1.
           S = F(M)
           S = SIGN(SQRT(ABS(S)),S)
           F(M) = 1.
           D(M) = S*D(M)*T
           U(H) = R*U(H)*S
           H = M
           R = T
           V = T
           IF ( JL .EQ. 1 ) GOTO 80
           IF ( JR .NE. 1 ) GOTO 100
           V = S
80         DO 90 I = 1,MQ
90              Q(I,M) = Q(I,M)*V
100        V = S
           IF ( JR .EQ. 2 ) GOTO 110
           IF ( JL .NE. 2 ) GOTO 130
           V = T
110        DO 120 I = 1,MP
120             P(I,M) = P(I,M)*V
130   CONTINUE
      RETURN
      END