C     SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S)                         SVD00100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1  SVD00200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974SVD00300
C          SINGULAR VALUE DECOMPOSITION ALSO TREATING RIGHT SIDE VECTOR.SVD00400
C                                                                       SVD00500
C     THE ARRAY S OCCUPIES 3*N CELLS.                                   SVD00600
C     A OCCUPIES M*N CELLS                                              SVD00700
C     B OCCUPIES M*NB CELLS.                                            SVD00800
C                                                                       SVD00900
C     SPECIAL SINGULAR VALUE DECOMPOSITION SUBROUTINE.                  SVD01000
C     WE HAVE THE M X N MATRIX A AND THE SYSTEM A*X=B TO SOLVE.         SVD01100
C     EITHER M .GE. N  OR  M .LT. N IS PERMITTED.                       SVD01200
C                       THE SINGULAR VALUE DECOMPOSITION                SVD01300
C     A = U*S*V**(T) IS MADE IN SUCH A WAY THAT ONE GETS                SVD01400
C        (1) THE MATRIX V IN THE FIRST N ROWS AND COLUMNS OF A.         SVD01500
C        (2) THE DIAGONAL MATRIX OF ORDERED SINGULAR VALUES IN          SVD01600
C            THE FIRST N CELLS OF THE ARRAY S(IP), IP .GE. 3*N.         SVD01700
C        (3) THE MATRIX PRODUCT U**(T)*B=G GETS PLACED BACK IN B.       SVD01800
C        (4) THE USER MUST COMPLETE THE SOLUTION AND DO HIS OWN         SVD01900
C            SINGULAR VALUE ANALYSIS.                                   SVD02000
C     *******                                                           SVD02100
C     GIVE SPECIAL                                                      SVD02200
C     TREATMENT TO ROWS AND COLUMNS WHICH ARE ENTIRELY ZERO.  THIS      SVD02300
C     CAUSES CERTAIN ZERO SING. VALS. TO APPEAR AS EXACT ZEROS RATHER   SVD02400
C     THAN AS ABOUT ETA TIMES THE LARGEST SING. VAL.   IT SIMILARLY     SVD02500
C     CLEANS UP THE ASSOCIATED COLUMNS OF U AND V.                      SVD02600
C     METHOD..                                                          SVD02700
C      1. EXCHANGE COLS OF A TO PACK NONZERO COLS TO THE LEFT.          SVD02800
C         SET N = NO. OF NONZERO COLS.                                  SVD02900
C         USE LOCATIONS A(1,NN),A(1,NN-1),...,A(1,N+1) TO RECORD THE    SVD03000
C         COL PERMUTATIONS.                                             SVD03100
C      2. EXCHANGE ROWS OF A TO PACK NONZERO ROWS TO THE TOP.           SVD03200
C         QUIT PACKING IF FIND N NONZERO ROWS.  MAKE SAME ROW EXCHANGES SVD03300
C         IN B.  SET M SO THAT ALL NONZERO ROWS OF THE PERMUTED A       SVD03400
C         ARE IN FIRST M ROWS.  IF M .LE. N THEN ALL M ROWS ARE         SVD03500
C         NONZERO.  IF M .GT. N THEN      THE FIRST N ROWS ARE KNOWN    SVD03600
C         TO BE NONZERO,AND ROWS N+1 THRU M MAY BE ZERO OR NONZERO.     SVD03700
C      3. APPLY ORIGINAL ALGORITHM TO THE M BY N PROBLEM.               SVD03800
C      4. MOVE PERMUTATION RECORD FROM A(,) TO S(I),I=N+1,...,NN.       SVD03900
C      5. BUILD V UP FROM  N BY N  TO  NN BY NN  BY PLACING ONES ON     SVD04000
C         THE DIAGONAL AND ZEROS ELSEWHERE.  THIS IS ONLY PARTLY DONE   SVD04100
C         EXPLICITLY.  IT IS COMPLETED DURING STEP 6.                   SVD04200
C      6. EXCHANGE ROWS OF V TO COMPENSATE FOR COL EXCHANGES OF STEP 2. SVD04300
C      7. PLACE ZEROS IN  S(I),I=N+1,NN  TO REPRESENT ZERO SING VALS.   SVD04400
C                                                                       SVD04500
      SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S)                         SVD04600
      DIMENSION A(MDA,NN),B(MDB,NB),S(NN,3)                             SVD04700
      ZERO=0.                                                           SVD04800
      ONE=1.                                                            SVD04900
C                                                                       SVD05000
C                             BEGIN.. SPECIAL FOR ZERO ROWS AND COLS.   SVD05100
C                                                                       SVD05200
C                             PACK THE NONZERO COLS TO THE LEFT         SVD05300
C                                                                       SVD05400
      N=NN                                                              SVD05500
      IF (N.LE.0.OR.MM.LE.0) RETURN                                     SVD05600
      J=N                                                               SVD05700
   10 CONTINUE                                                          SVD05800
          DO 20 I=1,MM                                                  SVD05900
          IF (A(I,J)) 50,20,50                                          SVD06000
   20     CONTINUE                                                      SVD06100
C                                                                       SVD06200
C         COL J  IS ZERO. EXCHANGE IT WITH COL N.                       SVD06300
C                                                                       SVD06400
      IF (J.EQ.N) GO TO 40                                              SVD06500
          DO 30 I=1,MM                                                  SVD06600
   30     A(I,J)=A(I,N)                                                 SVD06700
   40 CONTINUE                                                          SVD06800
      A(1,N)=J                                                          SVD06900
      N=N-1                                                             SVD07000
   50 CONTINUE                                                          SVD07100
      J=J-1                                                             SVD07200
      IF (J.GE.1) GO TO 10                                              SVD07300
C                             IF N=0 THEN A IS ENTIRELY ZERO AND SVD    SVD07400
C                             COMPUTATION CAN BE SKIPPED                SVD07500
      NS=0                                                              SVD07600
      IF (N.EQ.0) GO TO 240                                             SVD07700
C                             PACK NONZERO ROWS TO THE TOP              SVD07800
C                             QUIT PACKING IF FIND N NONZERO ROWS       SVD07900
      I=1                                                               SVD08000
      M=MM                                                              SVD08100
   60 IF (I.GT.N.OR.I.GE.M) GO TO 150                                   SVD08200
      IF (A(I,I)) 90,70,90                                              SVD08300
   70     DO 80 J=1,N                                                   SVD08400
          IF (A(I,J)) 90,80,90                                          SVD08500
   80     CONTINUE                                                      SVD08600
      GO TO 100                                                         SVD08700
   90 I=I+1                                                             SVD08800
      GO TO 60                                                          SVD08900
C                             ROW I IS ZERO                             SVD09000
C                             EXCHANGE ROWS I AND M                     SVD09100
  100 IF(NB.LE.0) GO TO 115                                             SVD09200
          DO 110 J=1,NB                                                 SVD09300
          T=B(I,J)                                                      SVD09400
          B(I,J)=B(M,J)                                                 SVD09500
  110     B(M,J)=T                                                      SVD09600
  115     DO 120 J=1,N                                                  SVD09700
  120     A(I,J)=A(M,J)                                                 SVD09800
      IF (M.GT.N) GO TO 140                                             SVD09900
          DO 130 J=1,N                                                  SVD10000
  130     A(M,J)=ZERO                                                   SVD10100
  140 CONTINUE                                                          SVD10200
C                             EXCHANGE IS FINISHED                      SVD10300
      M=M-1                                                             SVD10400
      GO TO 60                                                          SVD10500
C                                                                       SVD10600
  150 CONTINUE                                                          SVD10700
C                             END.. SPECIAL FOR ZERO ROWS AND COLUMNS   SVD10800
C                             BEGIN.. SVD ALGORITHM                     SVD10900
C     METHOD..                                                          SVD11000
C     (1)     REDUCE THE MATRIX TO UPPER BIDIAGONAL FORM WITH           SVD11100
C     HOUSEHOLDER TRANSFORMATIONS.                                      SVD11200
C          H(N)...H(1)AQ(1)...Q(N-2) = (D**T,0)**T                      SVD11300
C     WHERE D IS UPPER BIDIAGONAL.                                      SVD11400
C                                                                       SVD11500
C     (2)     APPLY H(N)...H(1) TO B.  HERE H(N)...H(1)*B REPLACES B    SVD11600
C     IN STORAGE.                                                       SVD11700
C                                                                       SVD11800
C     (3)     THE MATRIX PRODUCT W= Q(1)...Q(N-2) OVERWRITES THE FIRST  SVD11900
C     N ROWS OF A IN STORAGE.                                           SVD12000
C                                                                       SVD12100
C     (4)     AN SVD FOR D IS COMPUTED.  HERE K ROTATIONS RI AND PI ARE SVD12200
C     COMPUTED SO THAT                                                  SVD12300
C          RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SM)                SVD12400
C     TO WORKING ACCURACY.  THE SI ARE NONNEGATIVE AND NONINCREASING.   SVD12500
C     HERE RK...R1*B OVERWRITES B IN STORAGE WHILE                      SVD12600
C     A*P1**(T)...PK**(T)  OVERWRITES A IN STORAGE.                     SVD12700
C                                                                       SVD12800
C     (5)     IT FOLLOWS THAT,WITH THE PROPER DEFINITIONS,              SVD12900
C     U**(T)*B OVERWRITES B, WHILE V OVERWRITES THE FIRST N ROW AND     SVD13000
C     COLUMNS OF A.                                                     SVD13100
C                                                                       SVD13200
      L=MIN0(M,N)                                                       SVD13300
C             THE FOLLOWING LOOP REDUCES A TO UPPER BIDIAGONAL AND      SVD13400
C             ALSO APPLIES THE PREMULTIPLYING TRANSFORMATIONS TO B.     SVD13500
C                                                                       SVD13600
          DO 170 J=1,L                                                  SVD13700
          IF (J.GE.M) GO TO 160                                         SVD13800
          CALL H12 (1,J,J+1,M,A(1,J),1,T,A(1,J+1),1,MDA,N-J)            SVD13900
          CALL H12 (2,J,J+1,M,A(1,J),1,T,B,1,MDB,NB)                    SVD14000
  160     IF (J.GE.N-1) GO TO 170                                       SVD14100
          CALL H12 (1,J+1,J+2,N,A(J,1),MDA,S(J,3),A(J+1,1),MDA,1,M-J)   SVD14200
  170     CONTINUE                                                      SVD14300
C                                                                       SVD14400
C     COPY THE BIDIAGONAL MATRIX INTO THE ARRAY S() FOR QRBD.           SVD14500
C                                                                       SVD14600
      IF (N.EQ.1) GO TO 190                                             SVD14700
          DO 180 J=2,N                                                  SVD14800
          S(J,1)=A(J,J)                                                 SVD14900
  180     S(J,2)=A(J-1,J)                                               SVD15000
  190 S(1,1)=A(1,1)                                                     SVD15100
C                                                                       SVD15200
      NS=N                                                              SVD15300
      IF (M.GE.N) GO TO 200                                             SVD15400
      NS=M+1                                                            SVD15500
      S(NS,1)=ZERO                                                      SVD15600
      S(NS,2)=A(M,M+1)                                                  SVD15700
  200 CONTINUE                                                          SVD15800
C                                                                       SVD15900
C     CONSTRUCT THE EXPLICIT N BY N PRODUCT MATRIX, W=Q1*Q2*...*QL*I    SVD16000
C     IN THE ARRAY A().                                                 SVD16100
C                                                                       SVD16200
          DO 230 K=1,N                                                  SVD16300
          I=N+1-K                                                       SVD16400
          IF(I.GT.MIN0(M,N-2)) GO TO 210                                SVD16500
          CALL H12 (2,I+1,I+2,N,A(I,1),MDA,S(I,3),A(1,I+1),1,MDA,N-I)   SVD16600
  210         DO 220 J=1,N                                              SVD16700
  220         A(I,J)=ZERO                                               SVD16800
  230     A(I,I)=ONE                                                    SVD16900
C                                                                       SVD17000
C          COMPUTE THE SVD OF THE BIDIAGONAL MATRIX                     SVD17100
C                                                                       SVD17200
      CALL QRBD (IPASS,S(1,1),S(1,2),NS,A,MDA,N,B,MDB,NB)               SVD17300
C                                                                       SVD17400
      GO TO (240,310), IPASS                                            SVD17500
  240 CONTINUE                                                          SVD17600
      IF (NS.GE.N) GO TO 260                                            SVD17700
      NSP1=NS+1                                                         SVD17800
          DO 250 J=NSP1,N                                               SVD17900
  250     S(J,1)=ZERO                                                   SVD18000
  260 CONTINUE                                                          SVD18100
      IF (N.EQ.NN) RETURN                                               SVD18200
      NP1=N+1                                                           SVD18300
C                                  MOVE RECORD OF PERMUTATIONS          SVD18400
C                                  AND STORE ZEROS                      SVD18500
          DO 280 J=NP1,NN                                               SVD18600
          S(J,1)=A(1,J)                                                 SVD18700
              DO 270 I=1,N                                              SVD18800
  270         A(I,J)=ZERO                                               SVD18900
  280     CONTINUE                                                      SVD19000
C                             PERMUTE ROWS AND SET ZERO SINGULAR VALUES.SVD19100
          DO 300 K=NP1,NN                                               SVD19200
          I=S(K,1)                                                      SVD19300
          S(K,1)=ZERO                                                   SVD19400
              DO 290 J=1,NN                                             SVD19500
              A(K,J)=A(I,J)                                             SVD19600
  290         A(I,J)=ZERO                                               SVD19700
          A(I,K)=ONE                                                    SVD19800
  300     CONTINUE                                                      SVD19900
C                             END.. SPECIAL FOR ZERO ROWS AND COLUMNS   SVD20000
      RETURN                                                            SVD20100
  310 WRITE (6,320)                                                     SVD20200
      STOP                                                              SVD20300
  320 FORMAT (49H CONVERGENCE FAILURE IN QR BIDIAGONAL SVD ROUTINE)     SVD20400
      END                                                               SVD20500
C     SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)                QBD00100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 QBD00200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974QBD00300
C          QR ALGORITHM FOR SINGULAR VALUES OF A BIDIAGONAL MATRIX.     QBD00400
C                                                                       QBD00500
C     THE BIDIAGONAL MATRIX                                             QBD00600
C                                                                       QBD00700
C                       (Q1,E2,0...    )                                QBD00800
C                       (   Q2,E3,0... )                                QBD00900
C                D=     (       .      )                                QBD01000
C                       (         .   0)                                QBD01100
C                       (           .EN)                                QBD01200
C                       (          0,QN)                                QBD01300
C                                                                       QBD01400
C                 IS PRE AND POST MULTIPLIED BY                         QBD01500
C                 ELEMENTARY ROTATION MATRICES                          QBD01600
C                 RI AND PI SO THAT                                     QBD01700
C                                                                       QBD01800
C                 RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SN)         QBD01900
C                                                                       QBD02000
C                 TO WITHIN WORKING ACCURACY.                           QBD02100
C                                                                       QBD02200
C  1. EI AND QI OCCUPY E(I) AND Q(I) AS INPUT.                          QBD02300
C                                                                       QBD02400
C  2. RM...R1*C REPLACES 'C' IN STORAGE AS OUTPUT.                      QBD02500
C                                                                       QBD02600
C  3. V*P1**(T)...PM**(T) REPLACES 'V' IN STORAGE AS OUTPUT.            QBD02700
C                                                                       QBD02800
C  4. SI OCCUPIES Q(I) AS OUTPUT.                                       QBD02900
C                                                                       QBD03000
C  5. THE SI'S ARE NONINCREASING AND NONNEGATIVE.                       QBD03100
C                                                                       QBD03200
C     THIS CODE IS BASED ON THE PAPER AND 'ALGOL' CODE..                QBD03300
C REF..                                                                 QBD03400
C  1. REINSCH,C.H. AND GOLUB,G.H. 'SINGULAR VALUE DECOMPOSITION         QBD03500
C     AND LEAST SQUARES SOLUTIONS' (NUMER. MATH.), VOL. 14,(1970).      QBD03600
C                                                                       QBD03700
      SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)                QBD03800
      LOGICAL WNTV ,HAVERS,FAIL                                         QBD03900
      DIMENSION Q(NN),E(NN),V(MDV,NN),C(MDC,NCC)                        QBD04000
      ZERO=0.                                                           QBD04100
      ONE=1.                                                            QBD04200
      TWO=2.                                                            QBD04300
C                                                                       QBD04400
      N=NN                                                              QBD04500
      IPASS=1                                                           QBD04600
      IF (N.LE.0) RETURN                                                QBD04700
      N10=10*N                                                          QBD04800
      WNTV=NRV.GT.0                                                     QBD04900
      HAVERS=NCC.GT.0                                                   QBD05000
      FAIL=.FALSE.                                                      QBD05100
      NQRS=0                                                            QBD05200
      E(1)=ZERO                                                         QBD05300
      DNORM=ZERO                                                        QBD05400
           DO 10 J=1,N                                                  QBD05500
   10      DNORM=AMAX1(ABS(Q(J))+ABS(E(J)),DNORM)                       QBD05600
           DO 200 KK=1,N                                                QBD05700
           K=N+1-KK                                                     QBD05800
C                                                                       QBD05900
C     TEST FOR SPLITTING OR RANK DEFICIENCIES..                         QBD06000
C         FIRST MAKE TEST FOR LAST DIAGONAL TERM, Q(K), BEING SMALL.    QBD06100
   20       IF(K.EQ.1) GO TO 50                                         QBD06200
            IF(DIFF(DNORM+Q(K),DNORM)) 50,25,50                         QBD06300
C                                                                       QBD06400
C     SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO                QBD06500
C     TRANSFORM E(K) TO ZERO.                                           QBD06600
C                                                                       QBD06700
   25      CS=ZERO                                                      QBD06800
           SN=-ONE                                                      QBD06900
                DO 40 II=2,K                                            QBD07000
                I=K+1-II                                                QBD07100
                F=-SN*E(I+1)                                            QBD07200
                E(I+1)=CS*E(I+1)                                        QBD07300
                CALL G1 (Q(I),F,CS,SN,Q(I))                             QBD07400
C         TRANSFORMATION CONSTRUCTED TO ZERO POSITION (I,K).            QBD07500
C                                                                       QBD07600
                IF (.NOT.WNTV) GO TO 40                                 QBD07700
                     DO 30 J=1,NRV                                      QBD07800
   30                CALL G2 (CS,SN,V(J,I),V(J,K))                      QBD07900
C              ACCUMULATE RT. TRANSFORMATIONS IN V.                     QBD08000
C                                                                       QBD08100
   40           CONTINUE                                                QBD08200
C                                                                       QBD08300
C         THE MATRIX IS NOW BIDIAGONAL, AND OF LOWER ORDER              QBD08400
C         SINCE E(K) .EQ. ZERO..                                        QBD08500
C                                                                       QBD08600
   50           DO 60 LL=1,K                                            QBD08700
                L=K+1-LL                                                QBD08800
                IF(DIFF(DNORM+E(L),DNORM)) 55,100,55                    QBD08900
   55           IF(DIFF(DNORM+Q(L-1),DNORM)) 60,70,60                   QBD09000
   60           CONTINUE                                                QBD09100
C     THIS LOOP CAN'T COMPLETE SINCE E(1) = ZERO.                       QBD09200
C                                                                       QBD09300
           GO TO 100                                                    QBD09400
C                                                                       QBD09500
C         CANCELLATION OF E(L), L.GT.1.                                 QBD09600
   70      CS=ZERO                                                      QBD09700
           SN=-ONE                                                      QBD09800
                DO 90 I=L,K                                             QBD09900
                F=-SN*E(I)                                              QBD10000
                E(I)=CS*E(I)                                            QBD10100
                IF(DIFF(DNORM+F,DNORM)) 75,100,75                       QBD10200
   75           CALL G1 (Q(I),F,CS,SN,Q(I))                             QBD10300
                IF (.NOT.HAVERS) GO TO 90                               QBD10400
                     DO 80 J=1,NCC                                      QBD10500
   80                CALL G2 (CS,SN,C(I,J),C(L-1,J))                    QBD10600
   90           CONTINUE                                                QBD10700
C                                                                       QBD10800
C         TEST FOR CONVERGENCE..                                        QBD10900
  100      Z=Q(K)                                                       QBD11000
           IF (L.EQ.K) GO TO 170                                        QBD11100
C                                                                       QBD11200
C         SHIFT FROM BOTTOM 2 BY 2 MINOR OF B**(T)*B.                   QBD11300
           X=Q(L)                                                       QBD11400
           Y=Q(K-1)                                                     QBD11500
           G=E(K-1)                                                     QBD11600
           H=E(K)                                                       QBD11700
           F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)                        QBD11800
           G=SQRT(ONE+F**2)                                             QBD11900
           IF (F.LT.ZERO) GO TO 110                                     QBD12000
           T=F+G                                                        QBD12100
           GO TO 120                                                    QBD12200
  110      T=F-G                                                        QBD12300
  120      F=((X-Z)*(X+Z)+H*(Y/T-H))/X                                  QBD12400
C                                                                       QBD12500
C         NEXT QR SWEEP..                                               QBD12600
           CS=ONE                                                       QBD12700
           SN=ONE                                                       QBD12800
           LP1=L+1                                                      QBD12900
                DO 160 I=LP1,K                                          QBD13000
                G=E(I)                                                  QBD13100
                Y=Q(I)                                                  QBD13200
                H=SN*G                                                  QBD13300
                G=CS*G                                                  QBD13400
                CALL G1 (F,H,CS,SN,E(I-1))                              QBD13500
                F=X*CS+G*SN                                             QBD13600
                G=-X*SN+G*CS                                            QBD13700
                H=Y*SN                                                  QBD13800
                Y=Y*CS                                                  QBD13900
                IF (.NOT.WNTV) GO TO 140                                QBD14000
C                                                                       QBD14100
C              ACCUMULATE ROTATIONS (FROM THE RIGHT) IN 'V'             QBD14200
                     DO 130 J=1,NRV                                     QBD14300
  130                CALL G2 (CS,SN,V(J,I-1),V(J,I))                    QBD14400
  140           CALL G1 (F,H,CS,SN,Q(I-1))                              QBD14500
                F=CS*G+SN*Y                                             QBD14600
                X=-SN*G+CS*Y                                            QBD14700
                IF (.NOT.HAVERS) GO TO 160                              QBD14800
                     DO 150 J=1,NCC                                     QBD14900
  150                CALL G2 (CS,SN,C(I-1,J),C(I,J))                    QBD15000
C              APPLY ROTATIONS FROM THE LEFT TO                         QBD15100
C              RIGHT HAND SIDES IN 'C'..                                QBD15200
C                                                                       QBD15300
  160           CONTINUE                                                QBD15400
           E(L)=ZERO                                                    QBD15500
           E(K)=F                                                       QBD15600
           Q(K)=X                                                       QBD15700
           NQRS=NQRS+1                                                  QBD15800
           IF (NQRS.LE.N10) GO TO 20                                    QBD15900
C          RETURN TO 'TEST FOR SPLITTING'.                              QBD16000
C                                                                       QBD16100
           FAIL=.TRUE.                                                  QBD16200
C     ..                                                                QBD16300
C     CUTOFF FOR CONVERGENCE FAILURE. 'NQRS' WILL BE 2*N USUALLY.       QBD16400
  170      IF (Z.GE.ZERO) GO TO 190                                     QBD16500
           Q(K)=-Z                                                      QBD16600
           IF (.NOT.WNTV) GO TO 190                                     QBD16700
                DO 180 J=1,NRV                                          QBD16800
  180           V(J,K)=-V(J,K)                                          QBD16900
  190      CONTINUE                                                     QBD17000
C         CONVERGENCE. Q(K) IS MADE NONNEGATIVE..                       QBD17100
C                                                                       QBD17200
  200      CONTINUE                                                     QBD17300
      IF (N.EQ.1) RETURN                                                QBD17400
           DO 210 I=2,N                                                 QBD17500
           IF (Q(I).GT.Q(I-1)) GO TO 220                                QBD17600
  210      CONTINUE                                                     QBD17700
      IF (FAIL) IPASS=2                                                 QBD17800
      RETURN                                                            QBD17900
C     ..                                                                QBD18000
C     EVERY SINGULAR VALUE IS IN ORDER..                                QBD18100
  220      DO 270 I=2,N                                                 QBD18200
           T=Q(I-1)                                                     QBD18300
           K=I-1                                                        QBD18400
                DO 230 J=I,N                                            QBD18500
                IF (T.GE.Q(J)) GO TO 230                                QBD18600
                T=Q(J)                                                  QBD18700
                K=J                                                     QBD18800
  230           CONTINUE                                                QBD18900
           IF (K.EQ.I-1) GO TO 270                                      QBD19000
           Q(K)=Q(I-1)                                                  QBD19100
           Q(I-1)=T                                                     QBD19200
           IF (.NOT.HAVERS) GO TO 250                                   QBD19300
                DO 240 J=1,NCC                                          QBD19400
                T=C(I-1,J)                                              QBD19500
                C(I-1,J)=C(K,J)                                         QBD19600
  240           C(K,J)=T                                                QBD19700
  250      IF (.NOT.WNTV) GO TO 270                                     QBD19800
                DO 260 J=1,NRV                                          QBD19900
                T=V(J,I-1)                                              QBD20000
                V(J,I-1)=V(J,K)                                         QBD20100
  260           V(J,K)=T                                                QBD20200
  270      CONTINUE                                                     QBD20300
C         END OF ORDERING ALGORITHM.                                    QBD20400
C                                                                       QBD20500
      IF (FAIL) IPASS=2                                                 QBD20600
      RETURN                                                            QBD20700
      END                                                               QBD20800
