Power-Programmierung/CD1.mdf/fortran/library/ssp/sequence/dteas.for

Text File | 1985-11-29 | 4.9 KB | 207 lines

C
C     ..................................................................
C
C        SUBROUTINE DTEAS
C
C        PURPOSE
C           CALCULATE THE LIMIT OF A GIVEN SEQUENCE BY MEANS OF THE
C           EPSILON-ALGORITHM.
C
C        USAGE
C           CALL DTEAS(X,N,FIN,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X      - DOUBLE PRECISION VECTOR WHOSE COMPONENTS ARE TERMS
C                    OF THE GIVEN SEQUENCE. ON RETURN THE COMPONENTS OF
C                    VECTOR X ARE DESTROYED.
C           N      - DIMENSION OF INPUT VECTOR X.
C           FIN    - RESULTANT SCALAR IN DOUBLE PRECISION CONTAINING ON
C                    RETURN THE LIMIT OF THE GIVEN SEQUENCE.
C           EPS    - SINGLE PRECISION INPUT VALUE, WHICH SPECIFIES THE
C                    UPPER BOUND OF THE RELATIVE (ABSOLUTE) ERROR IF THE
C                    COMPONENTS OF X ARE ABSOLUTELY GREATER (LESS) THAN
C                    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 DTEAS(X,N,FIN,EPS,IER)
C
      DIMENSION X(1)
      DOUBLE PRECISION X,FIN,W1,W2,W3,W4,W5,W6,W7,T
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.D38
      W7=X(4)-X(3)
      IF(W7)3,4,3
    3 W1=1.D0/W7
C
    4 W5=1.D38
      W7=X(2)-X(1)
      IF(W7)5,6,5
    5 W5=1.D0/W7
C
    6 W4=X(3)-X(2)
      IF(W4)9,7,9
    7 W4=1.D38
      T=X(2)
      W2=X(3)
    8 W3=1.D38
      GO TO 17
C
    9 W4=1.D0/W4
C
      T=1.D38
      W7=W4-W5
      IF(W7)10,11,10
   10 T=X(2)+1.D0/W7
C
   11 W2=W1-W4
      IF(W2)15,12,15
   12 W2=1.D38
      IF(T-1.D38)13,14,14
   13 ISW2=1
   14 W3=W4
      GO TO 17
C
   15 W2=X(3)+1.D0/W2
      W7=W2-T
      IF(W7)16,8,16
   16 W3=W4+1.D0/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.D38
      W5=X(I-1)
      W7=X(I)-X(I-1)
      IF(W7)18,24,18
   18 W4=1.D0/W7
C
      IF(W1-1.D38)19,25,25
   19 W6=W4-W1
C
C        TEST FOR NECESSITY OF A SINGULAR RULE
C
      IF(DABS(W6)-DABS(W4)*1.D-12)20,20,22
   20 ISW2=1
      IF(W6)22,21,22
   21 W5=1.D38
      W6=W1
      IF(W2-1.D38)28,26,26
   22 W5=X(I-1)+1.D0/W6
C
C        FIRST TEST FOR LOSS OF SIGNIFICANCE
C
      IF(DABS(W5)-DABS(X(I-1))*1.D-10)23,24,24
   23 IF(W5)36,24,36
C
   24 W7=W5-W2
      IF(W7)27,25,27
   25 W6=1.D38
   26 ISW2=0
      X(IAUS)=W2
      GO TO 37
   27 W6=W1+1.D0/W7
   28 IF(ISW1-1)33,29,29
C
C        CALCULATE X(IAUS) WITH HELP OF SINGULAR RULE
C
   29 IF(W2-1.D38)30,32,32
   30 W7=W5/(W2-W5)+T/(W2-T)+X(I-2)/(X(I-2)-W2)
      IF(1.D0+W7)31,38,31
   31 X(IAUS)=W7*W2/(1.D0+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.D0/W7
C
C        SECOND TEST FOR LOSS OF SIGNIFICANCE
C
      IF(DABS(X(IAUS))-DABS(W2)*1.D-10)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.D38)39,38,38
   38 X(IAUS)=1.D38
      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
      HE1=DABS(X(I)-X(I+1))
      HE2=DABS(X(I+1))
      IF(HE1-EPS)44,44,42
   42 IF(HE2-1.)46,46,43
   43 IF(HE1-EPS*HE2)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