slaed4.f
来自「计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．」· F 代码 · 共 881 行 · 第 1/2 页
F
881 行
TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) END IF ELSE** (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)** We choose d(i+1) as origin.* ORGATI = .FALSE. DEL = D( IP1 ) - D( I ) A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) B = Z( IP1 )*Z( IP1 )*DEL IF( A.LT.ZERO ) THEN TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) ELSE TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) END IF END IF* IF( ORGATI ) THEN DO 130 J = 1, N DELTA( J ) = ( D( J )-D( I ) ) - TAU 130 CONTINUE ELSE DO 140 J = 1, N DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU 140 CONTINUE END IF IF( ORGATI ) THEN II = I ELSE II = I + 1 END IF IIM1 = II - 1 IIP1 = II + 1 OPS = OPS + 13*N + 6*( IIM1-IIP1 ) + 45** Evaluate PSI and the derivative DPSI* DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 150 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 150 CONTINUE ERRETM = ABS( ERRETM )** Evaluate PHI and the derivative DPHI* DPHI = ZERO PHI = ZERO DO 160 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 160 CONTINUE* W = RHOINV + PHI + PSI** W is the value of the secular function with* its ii-th element removed.* SWTCH3 = .FALSE. IF( ORGATI ) THEN IF( W.LT.ZERO ) $ SWTCH3 = .TRUE. ELSE IF( W.GT.ZERO ) $ SWTCH3 = .TRUE. END IF IF( II.EQ.1 .OR. II.EQ.N ) $ SWTCH3 = .FALSE.* TEMP = Z( II ) / DELTA( II ) DW = DPSI + DPHI + TEMP*TEMP TEMP = Z( II )*TEMP W = W + TEMP ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + $ THREE*ABS( TEMP ) + ABS( TAU )*DW** Test for convergence* IF( ABS( W ).LE.EPS*ERRETM ) THEN IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF** Calculate the new step* OPS = OPS + 14 NITER = NITER + 1 IF( .NOT.SWTCH3 ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )* $ ( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* $ ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - $ DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN OPS = OPS + 5 IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )* $ ( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )* $ ( DPSI+DPHI ) END IF END IF ETA = B / A ELSE IF( A.LE.ZERO ) THEN OPS = OPS + 8 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE OPS = OPS + 8 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF ELSE** Interpolation using THREE most relevant poles* OPS = OPS + 15 TEMP = RHOINV + PSI + PHI IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* $ ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* $ ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF ZZ( 2 ) = Z( II )*Z( II ) CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, $ INFO ) IF( INFO.NE.0 ) $ GO TO 250 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 + 18 + 7*N + 6*( IIM1-IIP1 ) IF( W*ETA.GE.ZERO ) THEN OPS = OPS + 1 ETA = -W / DW END IF TEMP = TAU + ETA DEL = ( D( IP1 )-D( I ) ) / TWO IF( ORGATI ) THEN IF( TEMP.GE.DEL ) THEN OPS = OPS + 1 ETA = DEL - TAU END IF IF( TEMP.LE.ZERO ) THEN OPS = OPS + 1 ETA = ETA / TWO END IF ELSE IF( TEMP.LE.-DEL ) THEN OPS = OPS + 1 ETA = -DEL - TAU END IF IF( TEMP.GE.ZERO ) THEN OPS = OPS + 1 ETA = ETA / TWO END IF END IF* PREW = W* 170 CONTINUE DO 180 J = 1, N DELTA( J ) = DELTA( J ) - ETA 180 CONTINUE** Evaluate PSI and the derivative DPSI* DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 190 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 190 CONTINUE ERRETM = ABS( ERRETM )** Evaluate PHI and the derivative DPHI* DPHI = ZERO PHI = ZERO DO 200 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 200 CONTINUE* TEMP = Z( 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+ETA )*DW* SWTCH = .FALSE. IF( ORGATI ) THEN IF( -W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. ELSE IF( W.GT.ABS( PREW ) / TEN ) $ SWTCH = .TRUE. END IF* TAU = TAU + ETA** Main loop to update the values of the array DELTA* ITER = NITER + 1* DO 240 NITER = ITER, MAXIT** Test for convergence* OPS = OPS + 1 IF( ABS( W ).LE.EPS*ERRETM ) THEN OPS = OPS + 1 IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF GO TO 250 END IF** Calculate the new step* IF( .NOT.SWTCH3 ) THEN OPS = OPS + 14 IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN C = W - DELTA( IP1 )*DW - $ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2 ELSE C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* $ ( Z( IP1 ) / DELTA( IP1 ) )**2 END IF ELSE TEMP = Z( II ) / DELTA( II ) IF( ORGATI ) THEN DPSI = DPSI + TEMP*TEMP ELSE DPHI = DPHI + TEMP*TEMP END IF C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI END IF A = ( DELTA( I )+DELTA( IP1 ) )*W - $ DELTA( I )*DELTA( IP1 )*DW B = DELTA( I )*DELTA( IP1 )*W IF( C.EQ.ZERO ) THEN IF( A.EQ.ZERO ) THEN OPS = OPS + 5 IF( .NOT.SWTCH ) THEN IF( ORGATI ) THEN A = Z( I )*Z( I ) + DELTA( IP1 )* $ DELTA( IP1 )*( DPSI+DPHI ) ELSE A = Z( IP1 )*Z( IP1 ) + $ DELTA( I )*DELTA( I )*( DPSI+DPHI ) END IF ELSE A = DELTA( I )*DELTA( I )*DPSI + $ DELTA( IP1 )*DELTA( IP1 )*DPHI END IF END IF OPS = OPS + 1 ETA = B / A ELSE IF( A.LE.ZERO ) THEN OPS = OPS + 8 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) ELSE OPS = OPS + 8 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) END IF ELSE** Interpolation using THREE most relevant poles* OPS = OPS + 2 TEMP = RHOINV + PSI + PHI IF( SWTCH ) THEN OPS = OPS + 10 C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI ELSE OPS = OPS + 14 IF( ORGATI ) THEN TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* $ ( ( DPSI-TEMP1 )+DPHI ) ELSE TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) TEMP1 = TEMP1*TEMP1 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* $ ( DPSI+( DPHI-TEMP1 ) ) ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) END IF END IF CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, $ INFO ) IF( INFO.NE.0 ) $ GO TO 250 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 + 7*N + 6*( IIM1-IIP1 ) + 18 IF( W*ETA.GE.ZERO ) THEN OPS = OPS + 1 ETA = -W / DW END IF TEMP = TAU + ETA DEL = ( D( IP1 )-D( I ) ) / TWO IF( ORGATI ) THEN IF( TEMP.GE.DEL ) THEN ETA = DEL - TAU OPS = OPS + 1 END IF IF( TEMP.LE.ZERO ) THEN ETA = ETA / TWO OPS = OPS + 1 END IF ELSE IF( TEMP.LE.-DEL ) THEN ETA = -DEL - TAU OPS = OPS + 1 END IF IF( TEMP.GE.ZERO ) THEN ETA = ETA / TWO OPS = OPS + 1 END IF END IF* DO 210 J = 1, N DELTA( J ) = DELTA( J ) - ETA 210 CONTINUE* TAU = TAU + ETA PREW = W** Evaluate PSI and the derivative DPSI* DPSI = ZERO PSI = ZERO ERRETM = ZERO DO 220 J = 1, IIM1 TEMP = Z( J ) / DELTA( J ) PSI = PSI + Z( J )*TEMP DPSI = DPSI + TEMP*TEMP ERRETM = ERRETM + PSI 220 CONTINUE ERRETM = ABS( ERRETM )** Evaluate PHI and the derivative DPHI* DPHI = ZERO PHI = ZERO DO 230 J = N, IIP1, -1 TEMP = Z( J ) / DELTA( J ) PHI = PHI + Z( J )*TEMP DPHI = DPHI + TEMP*TEMP ERRETM = ERRETM + PHI 230 CONTINUE* TEMP = Z( 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* 240 CONTINUE** Return with INFO = 1, NITER = MAXIT and not converged* INFO = 1 OPS = OPS + 1 IF( ORGATI ) THEN DLAM = D( I ) + TAU ELSE DLAM = D( IP1 ) + TAU END IF* END IF* 250 CONTINUE RETURN** End of SLAED4* END⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?