C     SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
C     TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE-HALL, 1974
C
C     MODIFIED AT SANDIA LABS, MAY 1977, TO INCREASE EFFICIENCY
C     THROUGH THE USE OF THE BASIC LINEAR ALGEBRA MODULE  SROT ,
C     WHICH PERFORMS THE GIVENS ROTATION APPLICATION STEP.
C
C          QR ALGORITHM FOR SINGULAR VALUES OF A BIDIAGONAL MATRIX.
C
C     THE BIDIAGONAL MATRIX
C
C                       (Q1,E2,0...    )
C                       (   Q2,E3,0... )
C                D=     (       .      )
C                       (         .   0)
C                       (           .EN)
C                       (          0,QN)
C
C                 IS PRE AND POST MULTIPLIED BY
C                 ELEMENTARY ROTATION MATRICES
C                 RI AND PI SO THAT
C
C                 RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SN)
C
C                 TO WITHIN WORKING ACCURACY.
C
C  1. EI AND QI OCCUPY E(I) AND Q(I) AS INPUT.
C
C  2. RM...R1*C REPLACES ^C^ IN STORAGE AS OUTPUT.
C
C  3. V*P1**(T)...PM**(T) REPLACES ^V^ IN STORAGE AS OUTPUT.
C
C  4. SI OCCUPIES Q(I) AS OUTPUT.
C
C  5. THE SI^S ARE NONINCREASING AND NONNEGATIVE.
C
C     THIS CODE IS BASED ON THE PAPER AND ^ALGOL^ CODE..
C REF..
C  1. REINSCH,C.H. AND GOLUB,G.H. ^SINGULAR VALUE DECOMPOSITION
C     AND LEAST SQUARES SOLUTIONS^ (NUMER. MATH.), VOL. 14,(1970).
C
      SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)
C
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
C
      LOGICAL WNTV ,HAVERS,FAIL
      DIMENSION Q(NN),E(NN),V(MDV,NN),C(MDC,NCC)
      ZERO=0.
      ONE=1.
      TWO=2.
C
      N=NN
      IPASS=1
      IF (N.LE.0) RETURN
      N10=10*N
      WNTV=NRV.GT.0
      HAVERS=NCC.GT.0
      FAIL=.FALSE.
      NQRS=0
      E(1)=ZERO
      DNORM=ZERO
           DO 10 J=1,N
   10      DNORM=MAX(ABS(Q(J))+ABS(E(J)),DNORM)
           DO 200 KK=1,N
           K=N+1-KK
C
C     TEST FOR SPLITTING OR RANK DEFICIENCIES..
C         FIRST MAKE TEST FOR LAST DIAGONAL TERM, Q(K), BEING SMALL.
   20       IF(K.EQ.1) GO TO 50
            IF(DIFF(DNORM+Q(K),DNORM)) 50,25,50
C
C     SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO
C     TRANSFORM E(K) TO ZERO.
C
   25      CS=ZERO
           SN=-ONE
                DO 40 II=2,K
                I=K+1-II
                F=-SN*E(I+1)
                E(I+1)=CS*E(I+1)
                CALL G1 (Q(I),F,CS,SN,Q(I))
C         TRANSFORMATION CONSTRUCTED TO ZERO POSITION (I,K).
C
                IF (.NOT.WNTV) GO TO 40
C     THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT.
C                    DO 30 J=1,NRV
C  30                CALL G2 (CS,SN,V(J,I),V(J,K))
                CALL SROT (NRV,V(1,I),1,V(1,K),1,CS,SN)
C              ACCUMULATE RT. TRANSFORMATIONS IN V.
C
   40           CONTINUE
C
C         THE MATRIX IS NOW BIDIAGONAL, AND OF LOWER ORDER
C         SINCE E(K) .EQ. ZERO..
C
   50           DO 60 LL=1,K
                L=K+1-LL
                IF(DIFF(DNORM+E(L),DNORM)) 55,100,55
   55           IF(DIFF(DNORM+Q(L-1),DNORM)) 60,70,60
   60           CONTINUE
C     THIS LOOP CAN^T COMPLETE SINCE E(1) = ZERO.
C
           GO TO 100
C
C         CANCELLATION OF E(L), L.GT.1.
   70      CS=ZERO
           SN=-ONE
                DO 90 I=L,K
                F=-SN*E(I)
                E(I)=CS*E(I)
                IF(DIFF(DNORM+F,DNORM)) 75,100,75
   75           CALL G1 (Q(I),F,CS,SN,Q(I))
                IF (.NOT.HAVERS) GO TO 90
C     THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT.
C                    DO 80 J=1,NCC
C  80                CALL G2 (CS,SN,C(I,J),C(L-1,J))
                CALL SROT (NCC,C(I,1),MDC,C(L-1,1),MDC,CS,SN)
   90           CONTINUE
C
C         TEST FOR CONVERGENCE..
  100      Z=Q(K)
           IF (L.EQ.K) GO TO 170
C
C         SHIFT FROM BOTTOM 2 BY 2 MINOR OF B**(T)*B.
           X=Q(L)
           Y=Q(K-1)
           G=E(K-1)
           H=E(K)
           F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)
           G=SQRT(ONE+F**2)
           IF (F.LT.ZERO) GO TO 110
           T=F+G
           GO TO 120
  110      T=F-G
  120      F=((X-Z)*(X+Z)+H*(Y/T-H))/X
C
C         NEXT QR SWEEP..
           CS=ONE
           SN=ONE
           LP1=L+1
                DO 160 I=LP1,K
                G=E(I)
                Y=Q(I)
                H=SN*G
                G=CS*G
                CALL G1 (F,H,CS,SN,E(I-1))
                F=X*CS+G*SN
                G=-X*SN+G*CS
                H=Y*SN
                Y=Y*CS
                IF (.NOT.WNTV) GO TO 140
C
C              ACCUMULATE ROTATIONS (FROM THE RIGHT) IN ^V^
C     THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT.
C                    DO 130 J=1,NRV
C 130                CALL G2 (CS,SN,V(J,I-1),V(J,I))
                CALL SROT (NRV,V(1,I-1),1,V(1,I),1,CS,SN)
  140           CALL G1 (F,H,CS,SN,Q(I-1))
                F=CS*G+SN*Y
                X=-SN*G+CS*Y
                IF (.NOT.HAVERS) GO TO 160
C     THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT.
C                    DO 150 J=1,NCC
C 150                CALL G2 (CS,SN,C(I-1,J),C(I,J))
                CALL SROT (NCC,C(I-1,1),MDC,C(I,1),MDC,CS,SN)
C              APPLY ROTATIONS FROM THE LEFT TO
C              RIGHT HAND SIDES IN ^C^..
C
  160           CONTINUE
           E(L)=ZERO
           E(K)=F
           Q(K)=X
           NQRS=NQRS+1
           IF (NQRS.LE.N10) GO TO 20
C          RETURN TO ^TEST FOR SPLITTING^.
C
           FAIL=.TRUE.
C     ..
C     CUTOFF FOR CONVERGENCE FAILURE. ^NQRS^ WILL BE 2*N USUALLY.
  170      IF (Z.GE.ZERO) GO TO 190
           Q(K)=-Z
           IF (.NOT.WNTV) GO TO 190
                DO 180 J=1,NRV
  180           V(J,K)=-V(J,K)
  190      CONTINUE
C         CONVERGENCE. Q(K) IS MADE NONNEGATIVE..
C
  200      CONTINUE
      IF (N.EQ.1) RETURN
           DO 210 I=2,N
           IF (Q(I).GT.Q(I-1)) GO TO 220
  210      CONTINUE
      IF (FAIL) IPASS=2
      RETURN
C     ..
C     EVERY SINGULAR VALUE IS IN ORDER..
  220      DO 270 I=2,N
           T=Q(I-1)
           K=I-1
                DO 230 J=I,N
                IF (T.GE.Q(J)) GO TO 230
                T=Q(J)
                K=J
  230           CONTINUE
           IF (K.EQ.I-1) GO TO 270
           Q(K)=Q(I-1)
           Q(I-1)=T
           IF (.NOT.HAVERS) GO TO 250
                DO 240 J=1,NCC
                T=C(I-1,J)
                C(I-1,J)=C(K,J)
  240           C(K,J)=T
  250      IF (.NOT.WNTV) GO TO 270
                DO 260 J=1,NRV
                T=V(J,I-1)
                V(J,I-1)=V(J,K)
  260           V(J,K)=T
  270      CONTINUE
C         END OF ORDERING ALGORITHM.
C
      IF (FAIL) IPASS=2
      RETURN
      END