home *** CD-ROM | disk | FTP | other *** search
- C
- C ..................................................................
- C
- C SUBROUTINE TEAS
- C
- C PURPOSE
- C CALCULATE THE LIMIT OF A GIVEN SEQUENCE BY MEANS OF THE
- C EPSILON-ALGORITHM.
- C
- C USAGE
- C CALL TEAS(X,N,FIN,EPS,IER)
- C
- C DESCRIPTION OF PARAMETERS
- C X - VECTOR WHOSE COMPONENTS ARE TERMS OF THE GIVEN
- C SEQUENCE. ON RETURN THE COMPONENTS OF VECTOR X
- C ARE DESTROYED.
- C N - DIMENSION OF INPUT VECTOR X.
- C FIN - RESULTANT SCALAR CONTAINING ON RETURN THE LIMIT
- C OF THE GIVEN SEQUENCE.
- C EPS - AN INPUT VALUE, WHICH SPECIFIES THE UPPER BOUND
- C OF THE RELATIVE (ABSOLUTE) ERROR IF THE COMPONENTS
- C OF X ARE ABSOLUTELY GREATER (LESS) THAN ONE.
- C CALCULATION IS TERMINATED AS SOON AS THREE TIMES IN
- C SUCCESSION THE RELATIVE (ABSOLUTE) DIFFERENCE
- C BETWEEN NEIGHBOURING TERMS IS NOT GREATER THAN EPS.
- C IER - RESULTANT ERROR PARAMETER CODED IN THE FOLLOWING
- C FORM
- C IER=0 - NO ERROR
- C IER=1 - REQUIRED ACCURACY NOT REACHED WITH
- C MAXIMAL NUMBER OF ITERATIONS
- C IER=-1 - INTEGER N IS LESS THAN TEN.
- C
- C REMARKS
- C NO ACTION BESIDES ERROR MESSAGE IN CASE N LESS THAN TEN.
- C THE CHARACTER OF THE GIVEN INFINITE SEQUENCE MUST BE
- C RECOGNIZABLE BY THOSE N COMPONENTS OF THE INPUT VECTOR X.
- C
- C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
- C NONE
- C
- C METHOD
- C THE CONVERGENCE OF THE GIVEN SEQUENCE IS ACCELERATED BY
- C MEANS OF THE E(2)-TRANSFORMATION, USED IN AN ITERATIVE WAY.
- C FOR REFERENCE, SEE
- C ALGORITHM 215,SHANKS, CACM 1963, NO. 11, PP. 662. AND
- C P. WYNN, SINGULAR RULES FOR CERTAIN NON-LINEAR ALGORITHMS
- C BIT VOL. 3, 1963, PP. 175-195.
- C
- C ..................................................................
- C
- SUBROUTINE TEAS(X,N,FIN,EPS,IER)
- C
- DIMENSION X(1)
- C
- C TEST ON WRONG INPUT PARAMETER N
- C
- NEW=N
- IF(NEW-10)1,2,2
- 1 IER=-1
- RETURN
- C
- C CALCULATE INITIAL VALUES FOR THE EPSILON ARRAY
- C
- 2 ISW1=0
- ISW2=0
- W1=1.E38
- W7=X(4)-X(3)
- IF(W7)3,4,3
- 3 W1=1./W7
- C
- 4 W5=1.E38
- W7=X(2)-X(1)
- IF(W7)5,6,5
- 5 W5=1./W7
- C
- 6 W4=X(3)-X(2)
- IF(W4)9,7,9
- 7 W4=1.E38
- T=X(2)
- W2=X(3)
- 8 W3=1.E38
- GO TO 17
- C
- 9 W4=1./W4
- C
- T=1.E38
- W7=W4-W5
- IF(W7)10,11,10
- 10 T=X(2)+1./W7
- C
- 11 W2=W1-W4
- IF(W2)15,12,15
- 12 W2=1.E38
- IF(T-1.E38)13,14,14
- 13 ISW2=1
- 14 W3=W4
- GO TO 17
- C
- 15 W2=X(3)+1./W2
- W7=W2-T
- IF(W7)16,8,16
- 16 W3=W4+1./W7
- C
- 17 ISW1=ISW2
- ISW2=0
- IMIN=4
- C
- C CALCULATE DIAGONALS OF THE EPSILON ARRAY IN A DO-LOOP
- C
- DO 40 I=5,NEW
- IAUS=I-IMIN
- W4=1.E38
- W5=X(I-1)
- W7=X(I)-X(I-1)
- IF(W7)18,24,18
- 18 W4=1./W7
- C
- IF(W1-1.E38)19,25,25
- 19 W6=W4-W1
- C
- C TEST FOR NECESSITY OF A SINGULAR RULE
- C
- IF(ABS(W6)-ABS(W4)*1.E-4)20,20,22
- 20 ISW2=1
- IF(W6)22,21,22
- 21 W5=1.E38
- W6=W1
- IF(W2-1.E38)28,26,26
- 22 W5=X(I-1)+1./W6
- C
- C FIRST TEST FOR LOSS OF SIGNIFICANCE
- C
- IF(ABS(W5)-ABS(X(I-1))*1.E-5)23,24,24
- 23 IF(W5)36,24,36
- C
- 24 W7=W5-W2
- IF(W7)27,25,27
- 25 W6=1.E38
- 26 ISW2=0
- X(IAUS)=W2
- GO TO 37
- 27 W6=W1+1./W7
- 28 IF(ISW1-1)33,29,29
- C
- C CALCULATE X(IAUS) WITH HELP OF SINGULAR RULE
- C
- 29 IF(W2-1.E38)30,32,32
- 30 W7=W5/(W2-W5)+T/(W2-T)+X(I-2)/(X(I-2)-W2)
- IF(1.+W7)31,38,31
- 31 X(IAUS)=W7*W2/(1.+W7)
- GO TO 39
- C
- 32 X(IAUS)=W5+T-X(I-2)
- GO TO 39
- C
- 33 W7=W6-W3
- IF(W7)34,38,34
- 34 X(IAUS)=W2+1./W7
- C
- C SECOND TEST FOR LOSS OF SIGNIFICANCE
- C
- IF(ABS(X(IAUS))-ABS(W2)*1.E-5)35,37,37
- 35 IF(X(IAUS))36,37,36
- C
- 36 NEW=IAUS-1
- ISW2=0
- GO TO 41
- C
- 37 IF(W2-1.E38)39,38,38
- 38 X(IAUS)=1.E38
- IMIN=I
- C
- 39 W1=W4
- T=W2
- W2=W5
- W3=W6
- ISW1=ISW2
- 40 ISW2=0
- C
- NEW=NEW-IMIN
- C
- C TEST FOR ACCURACY
- C
- 41 IEND=NEW-1
- DO 47 I=1,IEND
- W1=ABS(X(I)-X(I+1))
- W2=ABS(X(I+1))
- IF(W1-EPS)44,44,42
- 42 IF(W2-1.)46,46,43
- 43 IF(W1-EPS*W2)44,44,46
- 44 ISW2=ISW2+1
- IF(3-ISW2)45,45,47
- 45 FIN=X(I)
- IER=0
- RETURN
- C
- 46 ISW2=0
- 47 CONTINUE
- C
- IF(NEW-6)48,2,2
- 48 FIN=X(NEW)
- IER=1
- RETURN
- END