c subroutines TRIDIB, TINVIT from DEISPK -- these were added on 5 MAY 90 c (for documented version see specsubs.for) SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5) INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM DOUBLE PRECISION D(N),E(N),E2(N),W(M),RV4(N),RV5(N) DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON INTEGER IND(M) IERR = 0 TAG = 0 XU = D(1) X0 = D(1) U = 0.0D0 DO 40 I = 1, N X1 = U U = 0.0D0 IF (I .NE. N) U = DABS(E(I+1)) XU = DMIN1(D(I)-(X1+U),XU) X0 = DMAX1(D(I)+(X1+U),X0) IF (I .EQ. 1) GO TO 20 TST1 = DABS(D(I)) + DABS(D(I-1)) TST2 = TST1 + DABS(E(I)) IF (TST2 .GT. TST1) GO TO 40 20 E2(I) = 0.0D0 40 CONTINUE X1 = N X1 = X1 * EPSLON(DMAX1(DABS(XU),DABS(X0))) XU = XU - X1 T1 = XU X0 = X0 + X1 T2 = X0 P = 1 Q = N M1 = M11 - 1 IF (M1 .EQ. 0) GO TO 75 ISTURM = 1 50 V = X1 X1 = XU + (X0 - XU) * 0.5D0 IF (X1 .EQ. V) GO TO 980 GO TO 320 60 IF (S - M1) 65, 73, 70 65 XU = X1 GO TO 50 70 X0 = X1 GO TO 50 73 XU = X1 T1 = X1 75 M22 = M1 + M IF (M22 .EQ. N) GO TO 90 X0 = T2 ISTURM = 2 GO TO 50 80 IF (S - M22) 65, 85, 70 85 T2 = X1 90 Q = 0 R = 0 100 IF (R .EQ. M) GO TO 1001 TAG = TAG + 1 P = Q + 1 XU = D(P) X0 = D(P) U = 0.0D0 DO 120 Q = P, N X1 = U U = 0.0D0 V = 0.0D0 IF (Q .EQ. N) GO TO 110 U = DABS(E(Q+1)) V = E2(Q+1) 110 XU = DMIN1(D(Q)-(X1+U),XU) X0 = DMAX1(D(Q)+(X1+U),X0) IF (V .EQ. 0.0D0) GO TO 140 120 CONTINUE 140 X1 = EPSLON(DMAX1(DABS(XU),DABS(X0))) IF (EPS1 .LE. 0.0D0) EPS1 = -X1 IF (P .NE. Q) GO TO 180 IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940 M1 = P M2 = P RV5(P) = D(P) GO TO 900 180 X1 = X1 * (Q - P + 1) LB = DMAX1(T1,XU-X1) UB = DMIN1(T2,X0+X1) X1 = LB ISTURM = 3 GO TO 320 200 M1 = S + 1 X1 = UB ISTURM = 4 GO TO 320 220 M2 = S IF (M1 .GT. M2) GO TO 940 X0 = UB ISTURM = 5 DO 240 I = M1, M2 RV5(I) = UB RV4(I) = LB 240 CONTINUE K = M2 250 XU = LB DO 260 II = M1, K I = M1 + K - II IF (XU .GE. RV4(I)) GO TO 260 XU = RV4(I) GO TO 280 260 CONTINUE 280 IF (X0 .GT. RV5(K)) X0 = RV5(K) 300 X1 = (XU + X0) * 0.5D0 IF ((X0 - XU) .LE. DABS(EPS1)) GO TO 420 TST1 = 2.0D0 * (DABS(XU) + DABS(X0)) TST2 = TST1 + (X0 - XU) IF (TST2 .EQ. TST1) GO TO 420 320 S = P - 1 U = 1.0D0 DO 340 I = P, Q IF (U .NE. 0.0D0) GO TO 325 V = DABS(E(I)) / EPSLON(1.0D0) IF (E2(I) .EQ. 0.0D0) V = 0.0D0 GO TO 330 325 V = E2(I) / U 330 U = D(I) - X1 - V IF (U .LT. 0.0D0) S = S + 1 340 CONTINUE GO TO (60,80,200,220,360), ISTURM 360 IF (S .GE. K) GO TO 400 XU = X1 IF (S .GE. M1) GO TO 380 RV4(M1) = X1 GO TO 300 380 RV4(S+1) = X1 IF (RV5(S) .GT. X1) RV5(S) = X1 GO TO 300 400 X0 = X1 GO TO 300 420 RV5(K) = X1 K = K - 1 IF (K .GE. M1) GO TO 250 900 S = R R = R + M2 - M1 + 1 J = 1 K = M1 DO 920 L = 1, R IF (J .GT. S) GO TO 910 IF (K .GT. M2) GO TO 940 IF (RV5(K) .GE. W(L)) GO TO 915 DO 905 II = J, S I = L + S - II W(I+1) = W(I) IND(I+1) = IND(I) 905 CONTINUE 910 W(L) = RV5(K) IND(L) = TAG K = K + 1 GO TO 920 915 J = J + 1 920 CONTINUE 940 IF (Q .LT. N) GO TO 100 GO TO 1001 980 IERR = 3 * N + ISTURM 1001 LB = T1 UB = T2 RETURN END SUBROUTINE TINVIT(NM,N,D,E,E2,M,W,IND,Z, X IERR,RV1,RV2,RV3,RV4,RV6) INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M), X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) DOUBLE PRECISION U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,EPSLON, X PYTHAG INTEGER IND(M) IERR = 0 IF (M .EQ. 0) GO TO 1001 TAG = 0 ORDER = 1.0D0 - E2(1) Q = 0 100 P = Q + 1 DO 120 Q = P, N IF (Q .EQ. N) GO TO 140 IF (E2(Q+1) .EQ. 0.0D0) GO TO 140 120 CONTINUE 140 TAG = TAG + 1 S = 0 DO 920 R = 1, M IF (IND(R) .NE. TAG) GO TO 920 ITS = 1 X1 = W(R) IF (S .NE. 0) GO TO 510 XU = 1.0D0 IF (P .NE. Q) GO TO 490 RV6(P) = 1.0D0 GO TO 870 490 NORM = DABS(D(P)) IP = P + 1 DO 500 I = IP, Q 500 NORM = DMAX1(NORM, DABS(D(I))+DABS(E(I))) EPS2 = 1.0D-3 * NORM EPS3 = EPSLON(NORM) UK = Q - P + 1 EPS4 = UK * EPS3 UK = EPS4 / DSQRT(UK) S = P 505 GROUP = 0 GO TO 520 510 IF (DABS(X1-X0) .GE. EPS2) GO TO 505 GROUP = GROUP + 1 IF (ORDER * (X1 - X0) .LE. 0.0D0) X1 = X0 + ORDER * EPS3 520 V = 0.0D0 DO 580 I = P, Q RV6(I) = UK IF (I .EQ. P) GO TO 560 IF (DABS(E(I)) .LT. DABS(U)) GO TO 540 XU = U / E(I) RV4(I) = XU RV1(I-1) = E(I) RV2(I-1) = D(I) - X1 RV3(I-1) = 0.0D0 IF (I .NE. Q) RV3(I-1) = E(I+1) U = V - XU * RV2(I-1) V = -XU * RV3(I-1) GO TO 580 540 XU = E(I) / U RV4(I) = XU RV1(I-1) = U RV2(I-1) = V RV3(I-1) = 0.0D0 560 U = D(I) - X1 - XU * V IF (I .NE. Q) V = E(I+1) 580 CONTINUE IF (U .EQ. 0.0D0) U = EPS3 RV1(Q) = U RV2(Q) = 0.0D0 RV3(Q) = 0.0D0 600 DO 620 II = P, Q I = P + Q - II RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I) V = U U = RV6(I) 620 CONTINUE IF (GROUP .EQ. 0) GO TO 700 J = R DO 680 JJ = 1, GROUP 630 J = J - 1 IF (IND(J) .NE. TAG) GO TO 630 XU = 0.0D0 DO 640 I = P, Q 640 XU = XU + RV6(I) * Z(I,J) DO 660 I = P, Q 660 RV6(I) = RV6(I) - XU * Z(I,J) 680 CONTINUE 700 NORM = 0.0D0 DO 720 I = P, Q 720 NORM = NORM + DABS(RV6(I)) IF (NORM .GE. 1.0D0) GO TO 840 IF (ITS .EQ. 5) GO TO 830 IF (NORM .NE. 0.0D0) GO TO 740 RV6(S) = EPS4 S = S + 1 IF (S .GT. Q) S = P GO TO 780 740 XU = EPS4 / NORM DO 760 I = P, Q 760 RV6(I) = RV6(I) * XU 780 DO 820 I = IP, Q U = RV6(I) IF (RV1(I-1) .NE. E(I)) GO TO 800 U = RV6(I-1) RV6(I-1) = RV6(I) 800 RV6(I) = U - RV4(I) * RV6(I-1) 820 CONTINUE ITS = ITS + 1 GO TO 600 830 IERR = -R XU = 0.0D0 GO TO 870 840 U = 0.0D0 DO 860 I = P, Q 860 U = PYTHAG(U,RV6(I)) XU = 1.0D0 / U 870 DO 880 I = 1, N 880 Z(I,R) = 0.0D0 DO 900 I = P, Q 900 Z(I,R) = RV6(I) * XU X0 = X1 920 CONTINUE IF (Q .LT. N) GO TO 100 1001 RETURN END DOUBLE PRECISION FUNCTION EPSLON (X) DOUBLE PRECISION X DOUBLE PRECISION A,B,C,EPS A = 4.0D0/3.0D0 10 B = A - 1.0D0 C = B + B + B EPS = DABS(C-1.0D0) IF (EPS .EQ. 0.0D0) GO TO 10 EPSLON = EPS*DABS(X) RETURN END DOUBLE PRECISION FUNCTION PYTHAG(A,B) DOUBLE PRECISION A,B DOUBLE PRECISION P,R,S,T,U P = DMAX1(DABS(A),DABS(B)) IF (P .EQ. 0.0D0) GO TO 20 R = (DMIN1(DABS(A),DABS(B))/P)**2 10 CONTINUE T = 4.0D0 + R IF (T .EQ. 4.0D0) GO TO 20 S = R/T U = 1.0D0 + 2.0D0*S P = U*P R = (S/U)**2 * R GO TO 10 20 PYTHAG = P RETURN END SUBROUTINE COSSIN(N,C,S,C0,S0,F,DT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C(1),S(1) DATA PI/3.14159265358979D0/ C(1)=C0 S(1)=S0 CS=DCOS(2.D0*PI*F*DT) SN=DSIN(2.D0*PI*F*DT) DO 100 I=2,N C(I)=C(I-1)*CS-S(I-1)*SN 100 S(I)=C(I-1)*SN+S(I-1)*CS RETURN END SUBROUTINE REGRE(SR,SI,NF,NWIN,AMU) DIMENSION SR(4100,1),SI(4100,1),AMU(4100,1) COMMON/TAPSUM/B(10) C B IS FFT OF SLEPIAN EIGENTAPERS AT ZERO FREQ C SR SI ARE THE EIGENSPECTRA C AMU CONTAINS LINE FREQUENCY ESTIMATES AND F-TEST PARAMETER SUM=0. DO I=1,NWIN SUM=SUM+B(I)*B(I) END DO DO I=1,NF AMU(I,1)=0. AMU(I,2)=0. DO J=1,NWIN AMU(I,1)=AMU(I,1)+SR(I,J)*B(J) AMU(I,2)=AMU(I,2)+SI(I,J)*B(J) END DO AMU(I,1)=AMU(I,1)/SUM AMU(I,2)=AMU(I,2)/SUM SUM2=0. DO J=1,NWIN SUMR=SR(I,J)-AMU(I,1)*B(J) SUMI=SI(I,J)-AMU(I,2)*B(J) SUM2=SUM2+SUMR*SUMR+SUMI*SUMI END DO AMU(I,3)=(NWIN-1)*(AMU(I,2)**2+AMU(I,1)**2)*SUM/SUM2 END DO RETURN END SUBROUTINE FINDPJ(FMIN,FMAX,NF,F,P,NHW,NWIN,ipt) DIMENSION F(1),P(4100,1) character*80 ipt PMIN=3.00 ! ?% confidence cc PMIN=4.75 ! 95% confidence PRINT 901 901 FORMAT(' PER (kyr) FREQ (cyc/Myr) AMPLITUDE PHASE . $FTEST') OPEN(9,FILE='[.agu]'//ipt//'.ftest',status='new') WRITE(9,101) PMIN 101 FORMAT(E12.4,' IS THE PEAK THRESHHOLD') WRITE(9,901) I1=1 I2=2 OM=F(1)+F(2) 5 IF(OM.GE.FMIN) GO TO 10 I1=I2 I2=I2+1 OM=F(1)+(I2-1)*F(2) IF(I2.GE.NF) GO TO 987 GO TO 5 10 ISW=1 15 IF(P(I2,3).GT.P(I1,3)) ISW=2 GO TO (20,25),ISW 20 I1=I2 I2=I2+1 OM=F(1)+(I2-1)*F(2) IF(I2.GT.NF) GO TO 987 IF(OM.GT.FMAX) GO TO 987 GO TO 15 25 I1=I2 I2=I2+1 OM=F(1)+(I2-1)*F(2) IF(I2.GT.NF) GO TO 987 IF(OM.GT.FMAX) GO TO 987 IF(P(I2,3).GT.P(I1,3)) GO TO 25 IF(P(I1,3).LT.PMIN) GO TO 30 OMO=OM-F(2) omperiod=1000./omo PRINT 900,omperiod,OMO,(P(I1,IJK),IJK=1,3) WRITE(9,900) omperiod,OMO,(P(I1,IJK),IJK=1,3) 30 ISW=1 GO TO 20 900 FORMAT(5G14.6) 987 CLOSE(9) END cc gone!! SUBROUTINE GOUT(NF,NFZ,T,AMU,YLAB,IPT)