📄 dlasd4.f
字号:
END IF ELSE TEMP1 = Z( IIP1 ) / DTIIP TEMP1 = TEMP1*TEMP1 C = ( TEMP - DTIIM*( DPSI+DPHI ) ) - $ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1 IF( DPHI.LT.TEMP1 ) THEN OPS = OPS + DBLE( 2 ) ZZ( 1 ) = DTIIM*DTIIM*DPSI ELSE OPS = OPS + DBLE( 4 ) ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) ) END IF ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF OPS = OPS + DBLE( 2 ) ZZ( 2 ) = Z( II )*Z( II ) DD( 1 ) = DTIIM DD( 2 ) = DELTA( II )*WORK( II ) DD( 3 ) = DTIIP CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO ) IF( INFO.NE.0 ) $ GO TO 240 END IF** Note, eta should be positive if w is negative, and* eta should be negative otherwise. However,* if for some reason caused by roundoff, eta*w > 0,* we simply use one Newton step instead. This way* will guarantee eta*w < 0.* OPS = OPS + DBLE( 1 ) IF( W*ETA.GE.ZERO ) THEN OPS = OPS + DBLE( 1 ) ETA = -W / DW END IF OPS = OPS + DBLE( 8 ) IF( ORGATI ) THEN TEMP1 = WORK( I )*DELTA( I ) TEMP = ETA - TEMP1 ELSE TEMP1 = WORK( IP1 )*DELTA( IP1 ) TEMP = ETA - TEMP1 END IF IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN OPS = OPS + DBLE( 2 ) IF( W.LT.ZERO ) THEN ETA = ( SG2UB-TAU ) / TWO ELSE ETA = ( SG2LB-TAU ) / TWO END IF END IF* TAU = TAU + ETA ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )* PREW = W* OPS = OPS + DBLE( 1 + 2*N ) SIGMA = SIGMA + ETA DO 170 J = 1, N WORK( J ) = WORK( J ) + ETA DELTA( J ) = DELTA( J ) - ETA 170 CONTINUE** Evaluate PSI and the derivative DPSI* DPSI = ZERO PSI = ZERO ERRETM = ZERO OPS = OPS + DBLE( 7*IIM1 ) DO 180 J = 1, IIM1 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 180 CONTINUE ERRETM = ABS( ERRETM )** Evaluate PHI and the derivative DPHI* DPHI = ZERO PHI = ZERO OPS = OPS + DBLE( 7*(N-IIM1+1) ) DO 190 J = N, IIP1, -1 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 190 CONTINUE* OPS = OPS + DBLE( 19 ) TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU )*DW* IF( W.LE.ZERO ) THEN SG2LB = MAX( SG2LB, TAU ) ELSE SG2UB = MIN( SG2UB, TAU ) END IF* SWTCH = .FALSE. IF( ORGATI ) THEN IF( -W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. ELSE IF( W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. END IF** Main loop to update the values of the array DELTA and WORK* ITER = NITER + 1* DO 230 NITER = ITER, MAXIT** Test for convergence* OPS = OPS + DBLE( 1 ) IF( ABS( W ).LE.EPS*ERRETM ) THEN GO TO 240 END IF** Calculate the new step* IF( .NOT.SWTCH3 ) THEN OPS = OPS + DBLE( 2 ) DTIPSQ = WORK( IP1 )*DELTA( IP1 ) DTISQ = WORK( I )*DELTA( I ) IF( .NOT.SWTCH ) THEN OPS = OPS + DBLE( 6 ) IF( ORGATI ) THEN C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2 ELSE C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2 END IF ELSE OPS = OPS + DBLE( 8 ) TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) IF( ORGATI ) THEN DPSI = DPSI + TEMP*TEMP ELSE DPHI = DPHI + TEMP*TEMP END IF C = W - DTISQ*DPSI - DTIPSQ*DPHI END IF OPS = OPS + DBLE( 7 ) A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW B = DTIPSQ*DTISQ*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN OPS = OPS + DBLE( 5 ) IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ* $ ( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + $ DTISQ*DTISQ*( DPSI+DPHI ) END IF ELSE A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI END IF END IF OPS = OPS + DBLE( 1 ) ETA = B / A ELSE IF( A.LE.ZERO ) THEN OPS = OPS + DBLE( 8 ) ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE OPS = OPS + DBLE( 8 ) ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF ELSE** Interpolation using THREE most relevant poles* OPS = OPS + DBLE( 4 ) DTIIM = WORK( IIM1 )*DELTA( IIM1 ) DTIIP = WORK( IIP1 )*DELTA( IIP1 ) TEMP = RHOINV + PSI + PHI IF( SWTCH ) THEN OPS = OPS + DBLE( 8 ) C = TEMP - DTIIM*DPSI - DTIIP*DPHI ZZ( 1 ) = DTIIM*DTIIM*DPSI ZZ( 3 ) = DTIIP*DTIIP*DPHI ELSE IF( ORGATI ) THEN OPS = OPS + DBLE( 11 ) TEMP1 = Z( IIM1 ) / DTIIM TEMP1 = TEMP1*TEMP1 TEMP2 = ( D( IIM1 )-D( IIP1 ) )* $ ( D( IIM1 )+D( IIP1 ) )*TEMP1 C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) IF( DPSI.LT.TEMP1 ) THEN OPS = OPS + DBLE( 2 ) ZZ( 3 ) = DTIIP*DTIIP*DPHI ELSE OPS = OPS + DBLE( 4 ) ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI ) END IF ELSE OPS = OPS + DBLE( 10 ) TEMP1 = Z( IIP1 ) / DTIIP TEMP1 = TEMP1*TEMP1 TEMP2 = ( D( IIP1 )-D( IIM1 ) )* $ ( D( IIM1 )+D( IIP1 ) )*TEMP1 C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2 IF( DPHI.LT.TEMP1 ) THEN OPS = OPS + DBLE( 2 ) ZZ( 1 ) = DTIIM*DTIIM*DPSI ELSE OPS = OPS + DBLE( 4 ) ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) ) END IF OPS = OPS + DBLE( 1 ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF END IF OPS = OPS + DBLE( 1 ) DD( 1 ) = DTIIM DD( 2 ) = DELTA( II )*WORK( II ) DD( 3 ) = DTIIP CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO ) IF( INFO.NE.0 ) $ GO TO 240 END IF** Note, eta should be positive if w is negative, and* eta should be negative otherwise. However,* if for some reason caused by roundoff, eta*w > 0,* we simply use one Newton step instead. This way* will guarantee eta*w < 0.* OPS = OPS + DBLE( 1 ) IF( W*ETA.GE.ZERO ) THEN OPS = OPS + DBLE( 1 ) ETA = -W / DW END IF OPS = OPS + DBLE( 2 ) IF( ORGATI ) THEN TEMP1 = WORK( I )*DELTA( I ) TEMP = ETA - TEMP1 ELSE TEMP1 = WORK( IP1 )*DELTA( IP1 ) TEMP = ETA - TEMP1 END IF IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN OPS = OPS + DBLE( 2 ) IF( W.LT.ZERO ) THEN ETA = ( SG2UB-TAU ) / TWO ELSE ETA = ( SG2LB-TAU ) / TWO END IF END IF* OPS = OPS + DBLE( 6 ) TAU = TAU + ETA ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )* OPS = OPS + DBLE( 1 + 2*N ) SIGMA = SIGMA + ETA DO 200 J = 1, N WORK( J ) = WORK( J ) + ETA DELTA( J ) = DELTA( J ) - ETA 200 CONTINUE* PREW = W** Evaluate PSI and the derivative DPSI* DPSI = ZERO PSI = ZERO ERRETM = ZERO OPS = OPS + DBLE( 7*IIM1 ) DO 210 J = 1, IIM1 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 210 CONTINUE ERRETM = ABS( ERRETM )** Evaluate PHI and the derivative DPHI* DPHI = ZERO PHI = ZERO OPS = OPS + DBLE( 7*( IIM1-N+1 ) ) DO 220 J = N, IIP1, -1 TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 220 CONTINUE* OPS = OPS + DBLE( 19 ) TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = RHOINV + PHI + PSI + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU )*DW IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) $ SWTCH = .NOT.SWTCH* IF( W.LE.ZERO ) THEN SG2LB = MAX( SG2LB, TAU ) ELSE SG2UB = MIN( SG2UB, TAU ) END IF* 230 CONTINUE** Return with INFO = 1, NITER = MAXIT and not converged* INFO = 1* END IF* 240 CONTINUE RETURN** End of DLASD4* END⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -