SUBROUTINE D7RPLY(UU, VV, NZ, P, QP, K, QK, SVK) C VARIABLE-SHIFT K-POLYNOMIAL ITERATION FOR A C QUADRATIC FACTOR CONVERGES ONLY IF THE ZEROS ARE C EQUIMODULAR OR NEARLY SO. C UU,VV - COEFFICIENTS OF STARTING QUADRATIC C NZ - NUMBER OF ZERO FOUND C COMMON /P77PLY/ SR, SI, U, 1 V, A, B, C, D, A1, A2, A3, A6, A7, E, F, G, 2 H, SZR, SZI, LZR, LZI, ETA, ARE, MRE, N, NN C INTEGER N, NN INTEGER NZ, TYPE, I, J C DOUBLE PRECISION ETA, ARE, MRE DOUBLE PRECISION MP, OMP, EE, RELSTP, T, ZM C DOUBLE PRECISION P(1), QP(1), K(1), 1 QK(1), SVK(1), SR, SI, U, V, A, B, C, D, 2 A1, A2, A3, A6, A7, E, F, G, H, SZR, SZI, 3 LZR, LZI C DOUBLE PRECISION UI, VI, UU, VV LOGICAL TRIED C NZ = 0 TRIED = .FALSE. U = UU V = VV J = 0 C MAIN LOOP 10 CALL D6RPLY(1.D0, U, V, SZR, SZI, LZR, LZI) C RETURN IF ROOTS OF THE QUADRATIC ARE REAL AND NOT C CLOSE TO MULTIPLE OR NEARLY EQUAL AND OF OPPOSITE C SIGN IF (DABS(DABS(SZR)-DABS(LZR)).GT..01D0* 1 DABS(LZR)) RETURN C EVALUATE POLYNOMIAL BY QUADRATIC SYNTHETIC DIVISION CALL D8RPLY(NN, U, V, P, QP, A, B) MP = DABS(A-SZR*B) + DABS(SZI*B) C COMPUTE A RIGOROUS BOUND ON THE ROUNDING ERROR IN C EVALUTING P ZM = SQRT(ABS(V)) EE = 2.*ABS(QP(1)) T = -SZR*B DO 20 I=2,N EE = EE*ZM + ABS(QP(I)) 20 CONTINUE EE = EE*ZM + ABS(A)+T EE = (5.*MRE+4.*ARE)*EE - (5.*MRE+2.*ARE)* 2 (ABS(A+T)+ABS(B)*ZM) + 3 2.*ARE*ABS(T) C ITERATION HAS CONVERGED SUFFICIENTLY IF THE C POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND IF (MP.GT.20.*EE) GO TO 30 NZ = 2 RETURN 30 J = J + 1 C STOP ITERATION AFTER 20 STEPS IF (J.GT.20) RETURN IF (J.LT.2) GO TO 50 IF (RELSTP.GT..01 .OR. MP.LT.OMP .OR. TRIED) 1 GO TO 50 C A CLUSTER APPEARS TO BE STALLING THE CONVERGENCE. C FIVE FIXED SHIFT STEPS ARE TAKEN WITH A U,V CLOSE C TO THE CLUSTER IF (RELSTP.LT.ETA) RELSTP = ETA RELSTP = SQRT(RELSTP) U = U - U*RELSTP V = V + V*RELSTP CALL D8RPLY(NN, U, V, P, QP, A, B) DO 40 I=1,5 CALL D2RPLY(TYPE, K, QK) CALL D5RPLY(TYPE, P, QP, K, QK) 40 CONTINUE TRIED = .TRUE. J = 0 50 OMP = MP C CALCULATE NEXT K POLYNOMIAL AND NEW U AND V CALL D2RPLY(TYPE, K, QK) CALL D5RPLY(TYPE, P, QP, K, QK) CALL D2RPLY(TYPE, K, QK) CALL D4RPLY(TYPE, UI, VI, P, K) C IF VI IS ZERO THE ITERATION IS NOT CONVERGING IF (VI.EQ.0.D0) RETURN RELSTP = DABS((VI-V)/VI) U = UI V = VI GO TO 10 END .