MINICERN/mathlib/gen/c/r7dp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:57  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       SUBROUTINE C309R7(A,B,IBEG,INUM,XX,EPS)
  12 C
  13 C*******************************************************************
  14 C
  15 C  RCF converts polynomial A to the corresponding continued
  16 C         fraction, in 'normal'  form with coefficients B
  17 C         by the 'P algorithm' of Patry & Gupta
  18 C
  19 C   A(z) = A1/z + A2/z**3 + A3/z**5 + ... + An/z**(2n-1)
  20 C
  21 C   B(z) = B1/z+ B2/z+ B3/z+ .../(z+ Bn/z)
  22 C
  23 C  data:
  24 C   A     vector A(k), k=1,INUM         input
  25 C   B     vector B(k), k=IBEG,INUM      output
  26 C   IBEG  order of first coef. calc.    input
  27 C   INUM  order of A, even or odd       input
  28 C   XX    auxiliary vector of length .ge. length of vector B
  29 C         caller provides space for A,B,XX
  30 C   Note that neither of the first two terms A(1) A(2) should be zero
  31 C          & the user can start the calculation with any value of
  32 C          IBEG provided the c.f. coefs have been already
  33 C          calculated up to INUM = IBEG-1
  34 C          & the method breaks down as soon as the absolute value
  35 C          of a c.f. coef. is less than EPS.    At the time of the
  36 C          break up  INUM has been replaced by minus times the number
  37 C          of this coefficient.
  38 C   algorithm: J. Patry & S. Gupta, EIR-Bericht 247, November 1973
  39 C              Eidg. Institut fur Reaktorforschung
  40 C              Wuerenlingen, Switzerland
  41 C   see also:  Haenggi, Roesel & Trautmann,
  42 C              J. Comput. Phys., v. 137, (1980) 242-258
  43 C   note:      restart procedure modified by I.J.Thompson
  44 C
  45 C*******************************************************************
  46 C
  47       IMPLICIT COMPLEX*16(A-H,O-Z)
  48       DIMENSION A(100),B(100),XX(2,100)
  49       LOGICAL EVEN
  50       DOUBLE PRECISION EPS
  51
  52       IBN=INUM
  53       IF(IBEG .GT. 4) GO TO 50
  54       IF(IBEG .EQ. 4) GO TO 20
  55       B(1)=A(1)
  56       IF(IBN .GE. 2) B(2)=-A(2)/A(1)
  57       IF(IBN .LT. 3) RETURN
  58       X0=A(3)/A(2)
  59       XX(2,1)=B(2)
  60       XX(1,1)=-X0
  61       XX(1,2)=0
  62       B(3)=-X0-B(2)
  63       X0=-B(3)*A(2)
  64       M=3
  65       MP12=2
  66       EVEN=.TRUE.
  67       IF(IBN .LE. 3) RETURN
  68    20 IF(ABS(B(3)) .LT. EPS*ABS(X0)) THEN
  69        INUM=-M
  70        RETURN
  71       END IF
  72       M=4
  73    30 X1=A(M)
  74       M2M1=MP12
  75       MP12=M2M1+1
  76       IF(EVEN) MP12=M2M1
  77       DO 40 K = 2,MP12
  78    40 X1=X1+A(M-K+1)*XX(1,K-1)
  79       B(M)=-X1/X0
  80       IF(M .GE. IBN) RETURN
  81    50 IF(ABS(B(M)) .LT. EPS*ABS(X0)) THEN
  82        INUM=-M
  83        RETURN
  84       END IF
  85       DO 60 K = M2M1,2,-1
  86    60 XX(2,K)=XX(1,K)+B(M)*XX(2,K-1)
  87       XX(2,1)=XX(1,1)+B(M)
  88       DO 70 K = 1,M2M1
  89       X0=XX(2,K)
  90       XX(2,K)=XX(1,K)
  91    70 XX(1,K)=X0
  92       X0=X1
  93       XX(1,M2M1+1)=0
  94       M=M+1
  95       EVEN=.NOT.EVEN
  96       GO TO 30
  97       END
  98 #endif