SUBROUTINE SOQDSV(N0,A,B,WORK,INFO) IMPLICIT NONE INTEGER N0, INFO, MAXITER, IINFO REAL A(*), B(*), WORK(*) INTEGER N, M, I, J, K INTEGER INDRV1, INDRV2, INDRV3, INDRV4, INDRV5 REAL TMP1, TMP2, TMP3, TMP4, TMP5 REAL TAU, TAU1, TAU2 REAL SIGMA, SIGMA1, SIGMA2, DESIG, T, DESIG0 REAL C1, S1, EPS, TOL REAL ONE, ZERO, HALF, TWO, HUNDRD, CONST PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0, TWO = 2.0E0) PARAMETER (HUNDRD = 100.0E0, CONST=0.75E0) * INTRINSIC REAL EXTERNAL SLAMCH REAL SLAMCH * IF (N0 .EQ. 1) THEN A(1)=ABS(A(1)) INFO = 0 RETURN ENDIF IF (N0 .EQ. 2) THEN CALL SLAS2(A(1),B(1),A(2),A(2),A(1)) INFO = 0 RETURN ENDIF * A(1:N0)=ABS(A(1:N0)) B(1:N0-1)=ABS(B(1:N0-1)) * INDRV1 = 0 INDRV2 = INDRV1+N0 INDRV3 = INDRV2+N0 INDRV4 = INDRV3+N0 INDRV5 = INDRV4+N0 * EPS = SLAMCH( 'Precision' ) TOL = HUNDRD*EPS * M = 1 N = N0 * IF (A(M) .LT. A(N)) THEN K=(N-M+1)/2+M-1 DO J=M,K TMP1=A(J) A(J)=A(N+M-J) A(N+M-J)=TMP1 TMP1=B(J) B(J)=B(N-1+M-J) B(N-1+M-J)=TMP1 ENDDO END IF * TMP1 = A(M) DO J = M, N-3 IF (B(J) .LE. TOL*TMP1) THEN B(J) = -ZERO WORK(INDRV2+J) = -ZERO M = J+1 TMP1 = A(J+1) ELSE TMP1 = A(J+1)*(TMP1/(TMP1+B(J))) ENDIF ENDDO * IF (B(N-2) .LE. TOL*TMP1) THEN B(N-2) = ZERO TMP1 = A(N-1) ELSE TMP1 = A(N-1)*(TMP1/(TMP1+B(N-2))) ENDIF * IF (B(N-1) .LE. TOL*TMP1) THEN B(N-1) = ZERO ENDIF * B(N) = ZERO WORK(INDRV2+N) = ZERO * 3000 SIGMA = -B(N) SIGMA1 = SQRT(EPS)*SIGMA SIGMA2 = TOL*SIGMA DESIG = -WORK(INDRV2+N) MAXITER = 30*(N-M+1) DO I = 1, MAXITER * 15 IF (N-M+1 .EQ. 1) THEN CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0) GO TO 700 ENDIF IF (N-M+1 .EQ. 2) THEN CALL SLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1)) CALL SLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0) CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0) GO TO 700 ENDIF * IF (B(N-2) .LE. SIGMA2) THEN CALL SLAS2(A(N-1), B(N-1), A(N), A(N), A(N-1)) CALL SLARTG7(SIGMA,DESIG,A(N-1),A(N-1),DESIG0) CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0) N = N - 2 GO TO 15 ENDIF IF (B(N-1) .LE. SIGMA2) THEN CALL SLARTG7(SIGMA,DESIG,A(N),A(N),DESIG0) N = N - 1 GO TO 15 ENDIF * IF (A(M) .LT. A(N)) THEN K=(N-M+1)/2+M-1 DO J=M,K TMP1=A(J) A(J)=A(N+M-J) A(N+M-J)=TMP1 TMP1=B(J) B(J)=B(N-1+M-J) B(N-1+M-J)=TMP1 ENDDO END IF * CALL SLAS2(A(N-1), B(N-1), A(N), TAU, TMP3) TAU = MIN(TAU,A(N)) IF (TAU .LE. SIGMA1) GO TO 350 TAU2 = MINVAL(A(M:N-1)) IF (TAU2 .LE. SIGMA1) GO TO 350 IF (TAU2 .LE. TAU) THEN TAU1 = TAU2 GO TO 140 ELSE TAU1 = TAU ENDIF * TMP4 = A(M)-TAU TMP5 = A(M)+TAU IF (TMP4 .LE. ZERO) THEN GO TO 140 ELSE TMP3 = SQRT(TMP4)*SQRT(TMP5) ENDIF DO J = M, N-2 C1 = ONE/SQRT(ONE+(B(J)/TMP3)**2) TMP4 = C1*A(J+1)-TAU TMP5 = C1*A(J+1)+TAU IF (TMP4 .LE. ZERO) THEN GO TO 140 ELSE TMP3 = SQRT(TMP4)*SQRT(TMP5) ENDIF ENDDO TMP4 = A(N)/SQRT(ONE+(B(N-1)/TMP3)**2)-TAU IF (TMP4 .LT. ZERO .AND. TAU+TMP4 .EQ. TAU) TMP4=ZERO IF (TMP4 .LT. ZERO) THEN TAU2 = MAX(TMP4+TAU,ZERO) IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU TAU = TAU2 ENDIF IF (TAU .GT. ZERO) GO TO 125 140 WORK(INDRV5+M:INDRV5+N)=ONE/A(M:N) WORK(INDRV1+N) = WORK(INDRV5+N) TMP3 = WORK(INDRV1+N) DO J = N-1,M,-1 WORK(INDRV3+J) = B(J)*WORK(INDRV5+J) WORK(INDRV1+J) = WORK(INDRV5+J)+ $ WORK(INDRV3+J)*WORK(INDRV1+J+1) TMP3 = MAX(TMP3,WORK(INDRV1+J)) ENDDO TMP3 = ONE/TMP3 WORK(INDRV1+M) = (WORK(INDRV1+M)*TMP3)*WORK(INDRV5+M) TMP1 = WORK(INDRV1+M) DO J = M+1,N WORK(INDRV4+J) = B(J-1)*WORK(INDRV5+J) WORK(INDRV1+J) = (WORK(INDRV1+J)*TMP3)*WORK(INDRV5+J)+ $ WORK(INDRV4+J)*WORK(INDRV1+J-1) TMP1 = MAX(TMP1,WORK(INDRV1+J)) ENDDO TAU = ZERO IF (TMP1 .GT. ZERO) THEN TAU = MAX(TAU,SQRT(TMP3)/SQRT(TMP1)) ENDIF TMP1 = ONE/TMP1 WORK(INDRV1+N) = WORK(INDRV1+N)*TMP1 WORK(INDRV2+N) = WORK(INDRV1+N)*WORK(INDRV5+N) TMP3 = WORK(INDRV2+N) DO J = N-1,M,-1 WORK(INDRV1+J) = WORK(INDRV1+J)*TMP1 WORK(INDRV2+J) = WORK(INDRV1+J)*WORK(INDRV5+J)+ $ WORK(INDRV3+J)*WORK(INDRV2+J+1) TMP3 = MAX(TMP3,WORK(INDRV2+J)) ENDDO TMP3 = ONE/TMP3 TMP2 = (WORK(INDRV2+M)*TMP3)*WORK(INDRV5+M) TMP1 = TMP2/WORK(INDRV1+M) DO J = M+1,N TMP2 = (WORK(INDRV2+J)*TMP3)*WORK(INDRV5+J) $ +WORK(INDRV4+J)*TMP2 TMP1 = MAX(TMP1,TMP2/WORK(INDRV1+J)) ENDDO IF (TMP1 .GT. ZERO) THEN TAU = MAX(TAU,SQRT(TMP3)/SQRT(TMP1)) ENDIF TAU = MIN(TAU,TAU1) IF (TAU .GT. ZERO) GO TO 125 TAU=A(N)-HALF*B(N-1) IF (TAU .LE. ZERO) GO TO 350 DO J=N-1,M+1,-1 TMP2=A(J)-HALF*(B(J)+B(J-1)) IF (TMP2 .LE. ZERO) GO TO 350 TAU=MIN(TAU,TMP2) ENDDO TMP2=A(M)-HALF*B(M) IF (TMP2 .LE. ZERO) GO TO 350 TAU=MIN(TAU,TMP2) * 125 IF (TAU .LE. ZERO) GO TO 350 * TMP4 = A(M)-TAU TMP5 = A(M)+TAU IF (TMP4 .LE. ZERO) THEN TAU2 = MAX(TMP4+TAU,ZERO) IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU TAU = TAU2 GO TO 125 ELSE TMP3 = SQRT(TMP4)*SQRT(TMP5) ENDIF DO J = M, N-2 CALL SLARTG(TMP3,B(J),C1,S1,WORK(INDRV1+J)) WORK(INDRV2+J) = S1*A(J+1) TMP4 = C1*A(J+1)-TAU TMP5 = C1*A(J+1)+TAU IF (TMP4 .LE. ZERO) THEN TAU2 = MAX(TMP4+TAU,ZERO) IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU TAU = TAU2 GO TO 125 ELSE TMP3 = SQRT(TMP4)*SQRT(TMP5) ENDIF ENDDO CALL SLARTG(TMP3,B(N-1),C1,S1,WORK(INDRV1+N-1)) WORK(INDRV2+N-1) = S1*A(N) TMP4 = C1*A(N)-TAU TMP5 = C1*A(N)+TAU IF (TMP4 .LT. ZERO .AND. TAU+TMP4 .EQ. TAU) TMP4=ZERO IF (TMP4 .LT. ZERO) THEN TAU2 = MAX(TMP4+TAU,ZERO) IF (TAU2 .EQ. TAU) TAU2 = CONST*TAU IF (TAU2 .EQ. TAU) TAU2 = HALF*TAU TAU = TAU2 GO TO 125 ELSE TMP3 = SQRT(TMP4)*SQRT(TMP5) ENDIF IF (TMP3 .GE. A(N)) THEN TMP2 = TAU/A(N) IF (C1 .GE. TMP2) THEN TMP3 = A(N)*SQRT(C1-TMP2)*SQRT(C1+TMP2) ENDIF ENDIF WORK(INDRV1+N) = TMP3 * CALL SLARTG7(SIGMA,DESIG,TAU,T,DESIG0) * SIGMA = T SIGMA1 = SQRT(EPS)*SIGMA SIGMA2 = TOL*SIGMA DESIG = DESIG0 * GO TO 400 * 350 TMP1 = A(M) IF (TMP1 .LE. SIGMA1) TMP1=ZERO DO J = M, N-1 CALL SLARTG(TMP1,B(J),C1,S1,WORK(INDRV1+J)) WORK(INDRV2+J) = S1*A(J+1) TMP1 = C1*A(J+1) IF (TMP1 .LE. SIGMA1) TMP1=ZERO ENDDO WORK(INDRV1+N) = TMP1 * 400 TMP1 = WORK(INDRV1+M) TMP2 = TMP1 DO J = M, N-1 CALL SLARTG(TMP1,WORK(INDRV2+J),C1,S1,A(J)) B(J) = S1*WORK(INDRV1+J+1) TMP1 = C1*WORK(INDRV1+J+1) TMP2 = MIN(TMP2,TMP1) ENDDO A(N) = TMP1 IF (TMP2 .GT. SIGMA1) GO TO 450 TMP1 = WORK(INDRV1+M) IF (TMP1 .LE. SIGMA1) TMP1=ZERO DO J = M, N-1 CALL SLARTG(TMP1,WORK(INDRV2+J),C1,S1,A(J)) B(J) = S1*WORK(INDRV1+J+1) TMP1 = C1*WORK(INDRV1+J+1) IF (TMP1 .LE. SIGMA1) TMP1=ZERO ENDDO A(N) = TMP1 * 450 DO J = M, N-3 IF (B(J) .LE. SIGMA2 .OR. $ WORK(INDRV2+J) .LE. TOL*A(J)) THEN B(J) = -SIGMA WORK(INDRV2+J) = -DESIG M = J+1 ENDIF ENDDO IF (WORK(INDRV2+N-2) .LE. TOL*A(N-2)) THEN B(N-2) = ZERO ENDIF IF (WORK(INDRV2+N-1) .LE. TOL*A(N-1)) THEN B(N-1) = ZERO ENDIF * 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 CALL SLASRT( 'D', N0, A, IINFO ) * INFO = 0 RETURN * END