MINICERN/mathlib/gen/c/zerox64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:51  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DZEROX(A0,B0,EPS,MAXF,F,MODE)
  12 #endif
  13 #if !defined(CERNLIB_DOUBLE)
  14       FUNCTION RZEROX(A0,B0,EPS,MAXF,F,MODE)
  15 #endif
  16 C     Based on
  17 C
  18 C        J.C.P. Bus and T.J. Dekker, Two Efficient Algorithms with
  19 C        Guaranteed Convergence for Finding a Zero of a Function,
  20 C        ACM Trans. Math. Software 1 (1975) 330-345.
  21 C
  22 C        (MODE = 1: Algorithm M;    MODE = 2: Algorithm R)
  23 #include "gen/imp64.inc"
  24       CHARACTER NAME*(*)
  25       CHARACTER*80 ERRTXT
  26 #if defined(CERNLIB_DOUBLE)
  27       PARAMETER (NAME = 'DZEROX')
  28 #endif
  29 #if !defined(CERNLIB_DOUBLE)
  30       PARAMETER (NAME = 'RZEROX')
  31 #endif
  32       LOGICAL LMT
  33
  34       DIMENSION IM1(2),IM2(2),LMT(2)
  35
  36       PARAMETER (Z1 = 1, HALF = Z1/2)
  37
  38       DATA IM1 /2,3/, IM2 /-1,3/
  39
  40       IF(MODE .NE. 1 .AND. MODE .NE. 2) THEN
  41        C=0
  42        WRITE(ERRTXT,101) MODE
  43        CALL MTLPRT(NAME,'C200.1',ERRTXT)
  44        GO TO 99
  45       ENDIF
  46       FA=F(B0)
  47       FB=F(A0)
  48       IF(FA*FB .GT. 0) THEN
  49        C=0
  50        WRITE(ERRTXT,102) A0,B0
  51        CALL MTLPRT(NAME,'C200.2',ERRTXT)
  52        GO TO 99
  53       ENDIF
  54       ATL=ABS(EPS)
  55       B=A0
  56       A=B0
  57       LMT(2)=.TRUE.
  58       MF=2
  59     1 C=A
  60       FC=FA
  61     2 IE=0
  62     3 IF(ABS(FC) .LT. ABS(FB)) THEN
  63        IF(C .NE. A) THEN
  64         D=A
  65         FD=FA
  66        END IF
  67        A=B
  68        B=C
  69        C=A
  70        FA=FB
  71        FB=FC
  72        FC=FA
  73       END IF
  74       TOL=ATL*(1+ABS(C))
  75       H=HALF*(C+B)
  76       HB=H-B
  77       IF(ABS(HB) .GT. TOL) THEN
  78        IF(IE .GT. IM1(MODE)) THEN
  79         W=HB
  80        ELSE
  81         TOL=TOL*SIGN(Z1,HB)
  82         P=(B-A)*FB
  83         LMT(1)=IE .LE. 1
  84         IF(LMT(MODE)) THEN
  85          Q=FA-FB
  86          LMT(2)=.FALSE.
  87         ELSE
  88          FDB=(FD-FB)/(D-B)
  89          FDA=(FD-FA)/(D-A)
  90          P=FDA*P
  91          Q=FDB*FA-FDA*FB
  92         END IF
  93         IF(P .LT. 0) THEN
  94          P=-P
  95          Q=-Q
  96         END IF
  97         IF(IE .EQ. IM2(MODE)) P=P+P
  98         IF(P .EQ. 0 .OR. P .LE. Q*TOL) THEN
  99          W=TOL
 100         ELSEIF(P .LT. HB*Q) THEN
 101          W=P/Q
 102         ELSE
 103          W=HB
 104         END IF
 105        END IF
 106        D=A
 107        A=B
 108        FD=FA
 109        FA=FB
 110        B=B+W
 111        MF=MF+1
 112        IF(MF .GT. MAXF) THEN
 113         CALL MTLPRT(NAME,'C200.3','TOO MANY FUNCTION CALLS')
 114         GO TO 99
 115        ENDIF
 116        FB=F(B)
 117        IF(FB .EQ. 0 .OR. SIGN(Z1,FC) .EQ. SIGN(Z1,FB)) GO TO 1
 118        IF(W .EQ. HB) GO TO 2
 119        IE=IE+1
 120        GO TO 3
 121       END IF
 122 #if defined(CERNLIB_DOUBLE)
 123       DZEROX=C
 124 #endif
 125 #if !defined(CERNLIB_DOUBLE)
 126       RZEROX=C
 127 #endif
 128    99 CONTINUE
 129       RETURN
 130   101 FORMAT('MODE = ',I3,' ILLEGAL')
 131   102 FORMAT('F(A) AND F(B) HAVE THE SAME SIGN, A = ',1P,D15.8,
 132      1       ', B = ',D15.8)
 133       END
 134 #if !defined(CERNLIB_DOUBLE)
 135       FUNCTION ZEROX(A0,B0,EPS,MAXF,F,MODE)
 136       ZEROX=RZEROX(A0,B0,EPS,MAXF,F,MODE)
 137       END
 138 #endif