SUBROUTINE SSMALLSIGMA(FLAG,METHOD,M,N,A,B,TAUA,TAUB,TAU,
     $     WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7, RTMIN, RTMAX)

      IMPLICIT NONE
      
      REAL A(*), B(*), WORK2(*)
      INTEGER N, M, J, METHOD, FLAG
      INTEGER INDRV3, INDRV4, INDRV5, INDRV6, INDRV7
      REAL TMP1, TMP2, TMP3, TMP4, TMP5
      REAL TAUA(4), TAUB, TAU
      REAL C1, S1, T, RTMIN, RTMAX
      
      REAL ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0, TWO = 2.0E0)
      PARAMETER (CONST = 0.75E0)
      EXTERNAL SFMA0
      REAL SFMA0

      METHOD = 1
      IF (TAUB .LE. TAU) GO TO 160
*         
      TAUB = TAU
*     
      TMP4 = A(M)-TAU
      TMP5 = A(M)+TAU
      IF (TMP4 .LE. ZERO) THEN
         GO TO 160
      ELSE
         TMP3 = SQRT(TMP4)*SQRT(TMP5)
      ENDIF
      DO J = M, N-2
         CALL SLARTG2(TMP3,B(J),C1,S1,T,RTMIN,RTMAX)
         TMP4 = SFMA0(C1,A(J+1),-TAU)
         TMP5 = SFMA0(C1,A(J+1),TAU)
         IF (TMP4 .LE. ZERO) THEN
            GO TO 160
         ELSE
            TMP3 = SQRT(TMP4)*SQRT(TMP5)
         ENDIF
      ENDDO
      CALL SLARTG2(TMP3,B(N-1),C1,S1,T,RTMIN,RTMAX)
      TMP4 = SFMA0(C1,A(N),-TAU)
      IF (TMP4 .LT. ZERO) THEN
         TMP2 = C1*A(N)
         IF (TMP2 .LT. TAU) THEN
            IF (FLAG .EQ. 0) THEN
               TAUA(METHOD+1) = TMP2
            ELSE
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TAUA(METHOD+1) = MAX(TMP2,TMP3)
            ENDIF
         ELSE
            GO TO 160
         ENDIF
      ELSE
         TAUA(METHOD+1) = TAU
      ENDIF
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      IF (METHOD .GE. 2) RETURN
*
 160  WORK2(INDRV7+M:INDRV7+N) = ONE/A(M:N) 

      WORK2(INDRV5+N) = WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      TMP3 = ONE/TMP3
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)*WORK2(INDRV7+J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)*TMP3)*WORK2(INDRV7+J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      TMP1 = ONE/TMP1
      WORK2(INDRV5+N) = WORK2(INDRV5+N)*TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)*WORK2(INDRV7+N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)*TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)*WORK2(INDRV7+J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      TMP3 = ONE/TMP3
      TMP2 = (WORK2(INDRV6+M)*TMP3)*WORK2(INDRV7+M)
      TMP1 = TMP2/WORK2(INDRV5+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV6+J)*TMP3)*WORK2(INDRV7+J)
     $        +WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV5+J))
      ENDDO
      IF (TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(SQRT(TMP3/TMP1),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      IF (METHOD .GE. 2) RETURN
      
 170  WORK2(INDRV5+N) = ONE/A(N)
      TMP3 = WORK2(INDRV5+N)
      DO J = N-1,M,-1
         WORK2(INDRV3+J) = B(J)/A(J)
         WORK2(INDRV5+J) = ONE/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV5+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV5+J))
      ENDDO
      WORK2(INDRV5+M) = (WORK2(INDRV5+M)/TMP3)/A(M)
      TMP1 = WORK2(INDRV5+M)
      DO J = M+1,N
         WORK2(INDRV4+J) = B(J-1)/A(J)
         WORK2(INDRV5+J) = (WORK2(INDRV5+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*WORK2(INDRV5+J-1)
         TMP1 = MAX(TMP1,WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF
      
      WORK2(INDRV5+N) = WORK2(INDRV5+N)/TMP1
      WORK2(INDRV6+N) = WORK2(INDRV5+N)/A(N)
      TMP3 = WORK2(INDRV6+N)
      DO J = N-1,M,-1
         WORK2(INDRV5+J) = WORK2(INDRV5+J)/TMP1
         WORK2(INDRV6+J) = WORK2(INDRV5+J)/A(J)+
     $        WORK2(INDRV3+J)*WORK2(INDRV6+J+1)
         TMP3 = MAX(TMP3,WORK2(INDRV6+J))
      ENDDO
      TMP2 = (WORK2(INDRV6+M)/TMP3)/A(M)
      TMP1 = TMP2/WORK2(INDRV5+M)
      DO J = M+1,N
         TMP2 = (WORK2(INDRV6+J)/TMP3)/A(J)+
     $        WORK2(INDRV4+J)*TMP2
         TMP1 = MAX(TMP1,TMP2/WORK2(INDRV5+J))
      ENDDO
      IF (TMP3 .GT. ZERO .AND. TMP1 .GT. ZERO) THEN
         TAUA(METHOD+1) = MIN(ONE/SQRT(TMP1)/SQRT(TMP3),TAUB)
         IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
            METHOD = METHOD+1
         ENDIF
      ENDIF

      IF (METHOD .GE. 2) RETURN

      TAUA(METHOD+1) = A(N)-HALF*B(N-1)
      IF (TAUA(METHOD+1) .LE. ZERO) RETURN
      DO J=N-1,M+1,-1
         TMP2 = A(J)-HALF*(B(J)+B(J-1))
         IF (TMP2 .LE. ZERO) RETURN
         TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      ENDDO
      TMP2=A(M)-HALF*B(M)
      IF (TMP2 .LE. ZERO) RETURN
      TAUA(METHOD+1)=MIN(TAUA(METHOD+1),TMP2)
      IF (TAUA(METHOD+1) .GT. TAUA(METHOD)) THEN
         METHOD = METHOD+1
      ENDIF
      
      RETURN
      END SUBROUTINE SSMALLSIGMA

      SUBROUTINE SOQD(M,N,A,B,WORK2,
     $     INDRV1,INDRV2,INDRV5,INDRV6,INDRV7,INDRV8,SIGMA,SIGMA1,
     $     EPS,RTMIN,RTMAX,FLAG)

      IMPLICIT NONE
      INTEGER M, N, J, FLAG
      INTEGER INDRV1,INDRV2,INDRV5,INDRV6,INDRV7,INDRV8
      REAL A(*), B(*), WORK2(*)
      REAL TMP1, C1, S1, T
      REAL EPS, RTMIN, RTMAX, SIGMA, SIGMA1
      REAL ONE, ZERO, HALF, TWO
      PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0, TWO = 2.0E0)
      
      FLAG = -1
      TMP1 = A(M)
      IF (TMP1 .LE. SIGMA1) THEN
         CALL SLARTG2(SIGMA,TMP1,C1,S1,T,RTMIN,RTMAX)
         FLAG = M-1
         WORK2(INDRV1+M) = C1
         WORK2(INDRV2+M) = -S1
         TMP1 = ZERO
      ELSE
         WORK2(INDRV1+M) = ONE
         WORK2(INDRV2+M) = ZERO
      ENDIF
      DO J = M, N-1
         IF (B(J) .LE. EPS*TMP1) B(J) = ZERO
         CALL SLARTG2(TMP1,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
         WORK2(INDRV7+J) = C1
         WORK2(INDRV8+J) = S1
         WORK2(INDRV6+J) = S1*A(J+1)
         TMP1 = C1*A(J+1)
         IF (TMP1 .LE. SIGMA1) THEN
            CALL SLARTG2(SIGMA,TMP1,C1,S1,T,RTMIN,RTMAX)
            IF (FLAG .EQ. -1) FLAG = J
            WORK2(INDRV1+J+1) = C1
            WORK2(INDRV2+J+1) = -S1
            TMP1 = ZERO
         ELSE
            WORK2(INDRV1+J+1) = ONE
            WORK2(INDRV2+J+1) = ZERO
         ENDIF
      ENDDO
      WORK2(INDRV5+N) = TMP1      
      
      END SUBROUTINE SOQD
      
      SUBROUTINE SOQDS(UPLO,N0,A,B,SU,LDSU,SV,LDSV,WORK,WORK2,INFO)
      
      IMPLICIT NONE
      
      CHARACTER UPLO
      INTEGER N0, INFO, LDSU, LDSV, ISUB, MAXITER, FLAG
      REAL A(*), B(*), WORK2(*)
      REAL SU(LDSU, *), SV(LDSV, *), WORK(N0,*) 
      INTEGER N, M, I, J, K, OLDM, OLDN
      INTEGER INDRV1, INDRV2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7,
     $     INDRV8, METHOD
      REAL TMP1, TMP2, TMP3, TMP4, TMP5
      REAL TAU, TAUA(4), TAUB
      REAL SIGMA, SIGMA1, SIGMA2, DESIG, T, DESIG0, S
      REAL C1, S1, C2, S2, SMIN, EPS, TOL, SAFMIN
      REAL DMIN1, DMIN2
      REAL SAFMAX, RTMIN, RTMAX
      
      REAL ONE, ZERO, HALF, TWO, CONST
      PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0, TWO = 2.0E0)
      PARAMETER (CONST = 0.75E0)
      REAL TEN
      PARAMETER (TEN = 10.0D0)
*     
      INTRINSIC REAL
      EXTERNAL SLAMCH
      REAL SLAMCH
      EXTERNAL SFMA0
      REAL SFMA0
      LOGICAL LSAME
      EXTERNAL LSAME
*
      TAUA(1) = ZERO
*     
      INDRV1 = 0
      INDRV2 = INDRV1+N0
      INDRV3 = INDRV2+N0
      INDRV4 = INDRV3+N0
      INDRV5 = INDRV4+N0
      INDRV6 = INDRV5+N0
      INDRV7 = INDRV6+N0
      INDRV8 = INDRV7+N0
*     
      EPS = SLAMCH( 'Precision' )
      TOL = TEN*EPS
      SAFMIN = SLAMCH( 'S' )
      RTMIN = SQRT( SAFMIN )
      SAFMAX = ONE/SAFMIN
      RTMAX = SQRT( SAFMAX/2.0E0 )
*     
      IF ( LSAME( UPLO, 'U' ) ) THEN
         
         DO I = 1, N0-1
            IF (A(I) .LT. ZERO) THEN
               A(I) = -A(I)
               B(I) = -B(I)
               CALL SSCAL(N0,-ONE,SU(1,I),1)
            ENDIF
            IF (B(I) .LT. ZERO) THEN
               B(I) = -B(I)
               A(I+1) = -A(I+1)
               CALL SSCAL(N0,-ONE,SV(1,I+1),1)
            ENDIF
         ENDDO
         IF (A(N0) .LT. ZERO) THEN
            A(N0) = -A(N0)
            CALL SSCAL(N0,-ONE,SU(1,N0),1)
         ENDIF
         
         K=N0/2
         DO J=1,K
            TMP1=A(J)
            A(J)=A(N0+1-J)
            A(N0+1-J)=TMP1
            TMP1=B(J)
            B(J)=B(N0-J)
            B(N0-J)=TMP1
            CALL SSWAP(N0,SU(1,J),1,SU(1,N0+1-J),1)
            CALL SSWAP(N0,SV(1,J),1,SV(1,N0+1-J),1)
         ENDDO
            
      ELSE
*     
         DO I = 1, N0-1
            IF (A(I) .LT. ZERO) THEN
               A(I) = -A(I)
               B(I) = -B(I)
               CALL SSCAL(N0,-ONE,SV(1,I),1)
            ENDIF
            IF (B(I) .LT. ZERO) THEN
               B(I) = -B(I)
               A(I+1) = -A(I+1)
               CALL SSCAL(N0,-ONE,SU(1,I+1),1)
            ENDIF
         ENDDO
         IF (A(N0) .LT. ZERO) THEN
            A(N0) = -A(N0)
            CALL SSCAL(N0,-ONE,SV(1,N0),1)
         ENDIF
*     
      ENDIF
*
      OLDM = -1
      OLDN = -1
*     
      DO I = 1, N0
         DO J = 1, N0
            WORK(J,I)=ZERO
         ENDDO
      ENDDO
*     
      M = 1
      N = N0
*
      B(N) = ZERO
      WORK2(INDRV6+N) = ZERO
*     
 3000 SIGMA = -B(N)
      SIGMA1 = SQRT(EPS)*SIGMA
      SIGMA2 = TOL*SIGMA
      DESIG = -WORK2(INDRV6+N)
         
      MAXITER = 30*(N-M+1)
      DO I = 1, MAXITER
*     
 15      IF (N-M+1 .EQ. 1) THEN
            
            CALL SLARTG2(A(N),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)
            
            CALL SROT2(N0,SU(1,N),1,WORK(1,N),1,C1,S1)
            
            CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            
            GO TO 700
            
         ENDIF

         IF (N-M+1 .EQ. 2) THEN
            CALL SLASV2(A(N-1), B(N-1), A(N), A(N), A(N-1), S1, C1, 
     $           S2, C2)
            CALL SROT2(N0,SV(1,N-1),1,SV(1,N),1,C2,S2)
            CALL SROT2(N0,SU(1,N-1),1,SU(1,N),1,C1,S1)
            CALL SROT2(N0,WORK(1,N-1),1,WORK(1,N),1,C2,S2)
            CALL SLARTG2(A(N-1),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)
            CALL SROT2(N0,SU(1,N-1),1,WORK(1,N-1),1,C1,S1)
            CALL SLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0,RTMIN)
            CALL SLARTG2(A(N),SIGMA,C1,S1,TMP1,RTMIN,RTMAX)
            CALL SROT2(N0,SU(1,N),1,WORK(1,N),1,C1,S1)
            CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0,RTMIN)
            GO TO 700
         ENDIF

         IF (M .GT. OLDN .OR. N .LT. OLDM) THEN
            IF ( A(M) .LT. A(N) ) THEN

               TMP1 = A(N)
               DO J = N-1, M, -1
                  IF (B(J) .LE. MAX(SIGMA2,EPS*TMP1)) THEN
                     B(J) = -SIGMA
                     WORK2(INDRV6+J) = -DESIG
                     M = J+1
                     GO TO 15
                  ELSE
                     CALL SLARTG2(TMP1,B(J),C1,S1,TMP2,RTMIN,RTMAX)
                     TMP1 = A(J)*C1
                  ENDIF
               ENDDO
               
               DO J=M,N
                  WORK2(INDRV5+N+M-J)=A(J)
               ENDDO
               DO J=M,N-1
                  WORK2(INDRV6+N-1+M-J)=B(J)
               ENDDO

               K=(N-M+1)/2+M-1
               DO J=M,K
                  CALL SSWAP(N0,SU(1,J),1,SU(1,N+M-J),1)
                  CALL SSWAP(N0,WORK(1,J),1,WORK(1,N+M-J),1)
                  CALL SSWAP(N0,SV(1,J),1,SV(1,N+M-J),1)
               ENDDO

               OLDM = M
               OLDN = N
               
               TMP1 = WORK2(INDRV5+M)
               DMIN2 = TMP1
               DO J = M, N-3
                  IF (WORK2(INDRV6+J) .LE. EPS*TMP1)
     $                 WORK2(INDRV6+J) = ZERO
                  CALL SLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),
     $                 RTMIN,RTMAX)
                  WORK2(INDRV3+J) = C1
                  WORK2(INDRV4+J) = S1
                  B(J) = S1*WORK2(INDRV5+J+1)
                  TMP1 = C1*WORK2(INDRV5+J+1)
                  DMIN2 = MIN(DMIN2,TMP1)
               ENDDO

               IF (WORK2(INDRV6+N-2) .LE. EPS*TMP1)
     $              WORK2(INDRV6+N-2) = ZERO
               CALL SLARTG2(TMP1,WORK2(INDRV6+N-2),C1,S1,A(N-2),
     $              RTMIN,RTMAX)
               WORK2(INDRV3+N-2) = C1
               WORK2(INDRV4+N-2) = S1
               B(N-2) = S1*WORK2(INDRV5+N-1)
               TMP1 = C1*WORK2(INDRV5+N-1)
               DMIN1 = MIN(DMIN2,TMP1)

               IF (WORK2(INDRV6+N-1) .LE. EPS*TMP1)
     $              WORK2(INDRV6+N-1) = ZERO
               CALL SLARTG2(TMP1,WORK2(INDRV6+N-1),C1,S1,A(N-1),
     $              RTMIN,RTMAX)
               WORK2(INDRV3+N-1) = C1
               WORK2(INDRV4+N-1) = S1
               B(N-1) = S1*WORK2(INDRV5+N)
               TMP1 = C1*WORK2(INDRV5+N)

               A(N) = TMP1
*
               DO J = M, N-1
                  CALL SROT2(N0, WORK( 1, J), 1, WORK( 1, J+1), 1,
     $                 WORK2( INDRV3+J ), WORK2( INDRV4+J ))
               ENDDO
*     
               DO J = M, N-1
                  CALL SROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $                 WORK2( INDRV3+J ), WORK2( INDRV4+J ))
               ENDDO
               
               GO TO 400
            ENDIF
         ENDIF
*
         TMP1 = A(M)
         DO J = M, N-1
            IF (B(J) .LE. MAX(SIGMA2,EPS*TMP1)) THEN
               B(J) = -SIGMA
               WORK2(INDRV6+J) = -DESIG
               M = J+1
               GO TO 15
            ELSE
               CALL SLARTG2(TMP1,B(J),C1,S1,TMP2,RTMIN,RTMAX)
               TMP1 = A(J+1)*C1
            ENDIF
         ENDDO
*
         CALL SLAS2(A(N-1), B(N-1), A(N), TMP2, TMP3)
         TAU = MIN(TMP2,A(N))
         IF (TAU .LE. SIGMA1) THEN
            CALL SOQD(M,N,A,B,WORK2,
     $           INDRV1,INDRV2,INDRV5,INDRV6,INDRV7,INDRV8,
     $           SIGMA,SIGMA1,EPS,RTMIN,RTMAX,FLAG)
            IF (FLAG .GE. 0) GO TO 360
         ENDIF

         IF (OLDM .EQ. M .AND. OLDN .EQ. N) THEN
            TAUB = DMIN1
         ELSE IF (OLDM .EQ. M .AND. OLDN-1 .EQ. N) THEN
            TAUB = DMIN2
         ELSE
            TAUB = MINVAL(A(M:N-1))
         ENDIF
         IF (TAUB .LE. SIGMA1) THEN
            CALL SOQD(M,N,A,B,WORK2,
     $           INDRV1,INDRV2,INDRV5,INDRV6,INDRV7,INDRV8,
     $           SIGMA,SIGMA1,EPS,RTMIN,RTMAX,FLAG)
            IF (FLAG .GE. 0) GO TO 360
         ENDIF
*
         CALL SSMALLSIGMA(1,METHOD,M,N,A,B,TAUA,TAUB,TAU,
     $        WORK2, INDRV3, INDRV4, INDRV5, INDRV6, INDRV7,
     $        RTMIN, RTMAX)
         TAU = TAUA(METHOD)
*     
 125     IF (TAU .EQ. ZERO) GO TO 350
*     
         IF (SIGMA .EQ. ZERO) THEN
*     
            T = TAU
            DESIG0 = ZERO
*
            TMP4 = A(M)-TAU
            TMP5 = A(M)+TAU
            IF (TMP4 .LE. ZERO) THEN
               TMP2 = A(M)
               IF (TMP2 .GE. TAU) TMP2 = ZERO
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TMP2 = MAX(TMP2,TMP3)
               IF (METHOD .GE. 2) THEN
                  IF (TAUA(METHOD-1) .GE. TMP2) THEN
                     METHOD = METHOD-1
                     TAU = TAUA(METHOD)
                  ELSE
                     TAU = TMP2
                  ENDIF
               ELSE
                  TAU = TMP2
               ENDIF
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            CALL SLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
            
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = -S1
            DO J = M, N-2
               CALL SLARTG2(TMP3,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
               WORK2(INDRV7+J) = C1
               WORK2(INDRV8+J) = S1
               WORK2(INDRV6+J) = S1*A(J+1)
               TMP4 = SFMA0(C1,A(J+1),-TAU)
               TMP5 = SFMA0(C1,A(J+1),TAU)
               IF (TMP4 .LE. ZERO) THEN
                  TMP2 = C1*A(J+1)
                  IF (TMP2 .GE. TAU) TMP2 = ZERO
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
               CALL SLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = -S1
            ENDDO
            CALL SLARTG2(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1),
     $           RTMIN,RTMAX)
            WORK2(INDRV7+N-1) = C1
            WORK2(INDRV8+N-1) = S1
            WORK2(INDRV6+N-1) = S1*A(N)
            TMP4 = SFMA0(C1,A(N),-TAU)
            TMP5 = SFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO) THEN
               TMP2 = C1*A(N)
               IF (TMP2 .LT. TAU) THEN
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ENDIF
               TMP3 = ZERO
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF

            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = MIN(A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2),A(N))
               ELSE
                  TMP3 = A(N)
               ENDIF
            ENDIF
            
            CALL SLARTG2(TMP3,TAU,C1,S1,S,RTMIN,RTMAX)
            WORK2(INDRV1+N) = C1
            WORK2(INDRV2+N) = -S1
            WORK2(INDRV5+N) = TMP3
*     
         ELSE
*     
            CALL SLARTG7(SIGMA,DESIG,TAU,T,DESIG0,RTMIN)
*
            TMP4 = A(M)-TAU
            TMP5 = A(M)+TAU
            IF (TMP4 .LE. ZERO) THEN
               TMP2 = A(M)
               IF (TMP2 .GE. TAU) TMP2 = ZERO
               TMP3 = CONST*TAU
               IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
               TMP2 = MAX(TMP2,TMP3)
               IF (METHOD .GE. 2) THEN
                  IF (TAUA(METHOD-1) .GE. TMP2) THEN
                     METHOD = METHOD-1
                     TAU = TAUA(METHOD)
                  ELSE
                     TAU = TMP2
                  ENDIF
               ELSE
                  TAU = TMP2
               ENDIF
               GO TO 125
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF
            
            CALL SLARTG6(A(M),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
            WORK2(INDRV1+M) = C1
            WORK2(INDRV2+M) = S1
            DO J = M, N-2
               CALL SLARTG2(TMP3,B(J),C1,S1,WORK2(INDRV5+J),RTMIN,RTMAX)
               WORK2(INDRV7+J) = C1
               WORK2(INDRV8+J) = S1
               WORK2(INDRV6+J) = S1*A(J+1)
               
               TMP4 = SFMA0(C1,A(J+1),-TAU)
               TMP5 = SFMA0(C1,A(J+1),TAU)
               IF (TMP4 .LE. ZERO) THEN
                  TMP2 = C1*A(J+1)
                  IF (TMP2 .GE. TAU) TMP2 = ZERO
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ELSE
                  TMP3 = SQRT(TMP4)*SQRT(TMP5)
               ENDIF
               
               CALL SLARTG6(C1*A(J+1),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
               WORK2(INDRV1+J+1) = C1
               WORK2(INDRV2+J+1) = S1
            ENDDO
            CALL SLARTG2(TMP3,B(N-1),C1,S1,WORK2(INDRV5+N-1),
     $           RTMIN,RTMAX)
            WORK2(INDRV7+N-1) = C1
            WORK2(INDRV8+N-1) = S1
            WORK2(INDRV6+N-1) = S1*A(N)
            
            TMP4 = SFMA0(C1,A(N),-TAU)
            TMP5 = SFMA0(C1,A(N),TAU)
            IF (TMP4 .LT. ZERO) THEN
               TMP2 = C1*A(N)
               IF (TMP2 .LT. TAU) THEN
                  TMP3 = CONST*TAU
                  IF (TMP3 .GE. TAU) TMP3 = HALF*TAU
                  TMP2 = MAX(TMP2,TMP3)
                  IF (METHOD .GE. 2) THEN
                     IF (TAUA(METHOD-1) .GE. TMP2) THEN
                        METHOD = METHOD-1
                        TAU = TAUA(METHOD)
                     ELSE
                        TAU = TMP2
                     ENDIF
                  ELSE
                     TAU = TMP2
                  ENDIF
                  GO TO 125
               ENDIF
               TMP3 = ZERO
            ELSE
               TMP3 = SQRT(TMP4)*SQRT(TMP5)
            ENDIF

            IF (TMP3 .GE. A(N)) THEN
               TMP2 = TAU/A(N)
               IF (C1 .GE. TMP2) THEN
                  TMP3 = MIN(A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2),A(N))
               ELSE
                  TMP3 = A(N)
               ENDIF
            ENDIF
            
            IF (TMP3 .EQ. ZERO) THEN
               CALL SLARTG2(SIGMA,C1*A(N),C1,S1,S,RTMIN,RTMAX)
               S1 = -S1
            ELSE
               CALL SLARTG6(C1*A(N),TMP3,SIGMA,T,C1,S1,rtmin,rtmax)
            ENDIF
            WORK2(INDRV1+N) = C1
            WORK2(INDRV2+N) = S1
            WORK2(INDRV5+N) = TMP3
*     
         ENDIF
*     
         SIGMA = T
         SIGMA1 = SQRT(EPS)*SIGMA
         SIGMA2 = TOL*SIGMA
         DESIG = DESIG0
*
         OLDM = M
         OLDN = N
*
         TMP1 = WORK2(INDRV5+M)
         DMIN2 = TMP1
         DO J = M, N-3
            IF (WORK2(INDRV6+J) .LE. EPS*TMP1) WORK2(INDRV6+J) = ZERO
            CALL SLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),RTMIN,RTMAX)
            WORK2(INDRV3+J) = C1
            WORK2(INDRV4+J) = S1
            B(J) = S1*WORK2(INDRV5+J+1)
            TMP1 = C1*WORK2(INDRV5+J+1)
            DMIN2 = MIN(DMIN2,TMP1)
         ENDDO

         IF (WORK2(INDRV6+N-2) .LE. EPS*TMP1) WORK2(INDRV6+N-2) = ZERO
         CALL SLARTG2(TMP1,WORK2(INDRV6+N-2),C1,S1,A(N-2),RTMIN,RTMAX)
         WORK2(INDRV3+N-2) = C1
         WORK2(INDRV4+N-2) = S1
         B(N-2) = S1*WORK2(INDRV5+N-1)
         TMP1 = C1*WORK2(INDRV5+N-1)
         DMIN1 = MIN(DMIN2,TMP1)

         IF (WORK2(INDRV6+N-1) .LE. EPS*TMP1) WORK2(INDRV6+N-1) = ZERO
         CALL SLARTG2(TMP1,WORK2(INDRV6+N-1),C1,S1,A(N-1),RTMIN,RTMAX)
         WORK2(INDRV3+N-1) = C1
         WORK2(INDRV4+N-1) = S1
         B(N-1) = S1*WORK2(INDRV5+N)
         TMP1 = C1*WORK2(INDRV5+N)
         
         A(N) = TMP1
*
         CALL SROT2(N0,SU(1,M),1,WORK(1,M),1,
     $        WORK2(INDRV1+M),WORK2(INDRV2+M))
         DO J = M, N-1
            CALL SROT2(N0,SU(1,J),1,SU(1,J+1),1,WORK2(INDRV7+J),
     $           WORK2(INDRV8+J))
            CALL SROT2(N0,SU(1,J+1),1,WORK(1,J+1),1,
     $           WORK2(INDRV1+J+1),WORK2(INDRV2+J+1))
            CALL SROT2(N0,WORK(1,J),1,WORK(1,J+1),1,
     $           WORK2(INDRV3+J),WORK2(INDRV4+J))
         ENDDO
*
         DO J = M, N-1
            CALL SROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*     
         GO TO 400
*     
 350     CALL SOQD(M,N,A,B,WORK2,
     $        INDRV1,INDRV2,INDRV5,INDRV6,INDRV7,INDRV8,
     $        SIGMA,SIGMA1,EPS,RTMIN,RTMAX,FLAG)
*
 360     OLDM = M
         OLDN = N
*     
         TMP1 = WORK2(INDRV5+M)
         DMIN2 = TMP1
         DO J = M, N-3
            IF (WORK2(INDRV6+J) .LE. EPS*TMP1) WORK2(INDRV6+J) = ZERO
            CALL SLARTG2(TMP1,WORK2(INDRV6+J),C1,S1,A(J),RTMIN,RTMAX)
            WORK2(INDRV3+J) = C1
            WORK2(INDRV4+J) = S1
            B(J) = S1*WORK2(INDRV5+J+1)
            TMP1 = C1*WORK2(INDRV5+J+1)
            DMIN2 = MIN(DMIN2,TMP1)
         ENDDO

         IF (WORK2(INDRV6+N-2) .LE. EPS*TMP1) WORK2(INDRV6+N-2) = ZERO
         CALL SLARTG2(TMP1,WORK2(INDRV6+N-2),C1,S1,A(N-2),RTMIN,RTMAX)
         WORK2(INDRV3+N-2) = C1
         WORK2(INDRV4+N-2) = S1
         B(N-2) = S1*WORK2(INDRV5+N-1)
         TMP1 = C1*WORK2(INDRV5+N-1)
         DMIN1 = MIN(DMIN2,TMP1)

         IF (WORK2(INDRV6+N-1) .LE. EPS*TMP1) WORK2(INDRV6+N-1) = ZERO
         CALL SLARTG2(TMP1,WORK2(INDRV6+N-1),C1,S1,A(N-1),RTMIN,RTMAX)
         WORK2(INDRV3+N-1) = C1
         WORK2(INDRV4+N-1) = S1
         B(N-1) = S1*WORK2(INDRV5+N)
         TMP1 = C1*WORK2(INDRV5+N)

         A(N) = TMP1
*     
         IF (FLAG .EQ. -1) THEN

            DO J = M, N-1
               CALL SROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

         ELSE

            DO J = M, FLAG
               CALL SROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

            CALL SROT2(N0,SU(1,FLAG+1),1,WORK(1,FLAG+1),1,
     $           WORK2(INDRV1+FLAG+1),WORK2(INDRV2+FLAG+1))

            DO J = FLAG+1, N-1
               CALL SROT2(N0, SU( 1, J), 1, SU( 1, J+1), 1,
     $              WORK2( INDRV7+J ), WORK2( INDRV8+J ))
            ENDDO

         ENDIF
*
         DO J = M, N-1
            CALL SROT2(N0, WORK( 1, J), 1, WORK( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*     
         DO J = M, N-1
            CALL SROT2(N0, SV( 1, J), 1, SV( 1, J+1), 1,
     $           WORK2( INDRV3+J ), WORK2( INDRV4+J ))
         ENDDO
*     
 400  ENDDO
      INFO = 2
      RETURN
      
 700  N = M-1
      IF (N .LE. 0) THEN
         GO TO 4000
      ENDIF
*     
      DO J = M-2,1,-1
         IF (B(J) .LE. ZERO) THEN
            M = J+1
            GO TO 3000
         ENDIF
      ENDDO
      M = 1
      GO TO 3000
*     
 4000 DO I = 1, N0 - 1
*     
*     Scan for smallest D(I)
*     
         ISUB = 1
         SMIN = A( 1 )
         DO J = 2, N0 + 1 - I
            IF( A( J ).LE.SMIN ) THEN
               ISUB = J
               SMIN = A( J )
            END IF
         ENDDO
         IF( ISUB.NE.N0+1-I ) THEN
*     
*     Swap singular values and vectors
*     
            A( ISUB ) = A( N0+1-I )
            A( N0+1-I ) = SMIN
            CALL SSWAP( N0, SV( 1, ISUB ), 1, SV( 1, N0+1-I ), 1 ) 
            CALL SSWAP( N0, SU( 1, ISUB ), 1, SU( 1, N0+1-I ), 1 )
         END IF
      ENDDO
*
 30   INFO = 0
      RETURN
*     
      END