[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / r7dp.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:01:57  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if defined(CERNLIB_DOUBLE)
      SUBROUTINE C309R7(A,B,IBEG,INUM,XX,EPS)
C
C*******************************************************************
C
C  RCF converts polynomial A to the corresponding continued
C         fraction, in 'normal'  form with coefficients B
C         by the 'P algorithm' of Patry & Gupta
C
C   A(z) = A1/z + A2/z**3 + A3/z**5 + ... + An/z**(2n-1)
C
C   B(z) = B1/z+ B2/z+ B3/z+ .../(z+ Bn/z)
C
C  data:
C   A     vector A(k), k=1,INUM         input
C   B     vector B(k), k=IBEG,INUM      output
C   IBEG  order of first coef. calc.    input
C   INUM  order of A, even or odd       input
C   XX    auxiliary vector of length .ge. length of vector B
C         caller provides space for A,B,XX
C   Note that neither of the first two terms A(1) A(2) should be zero
C          & the user can start the calculation with any value of
C          IBEG provided the c.f. coefs have been already
C          calculated up to INUM = IBEG-1
C          & the method breaks down as soon as the absolute value
C          of a c.f. coef. is less than EPS.    At the time of the
C          break up  INUM has been replaced by minus times the number
C          of this coefficient.
C   algorithm: J. Patry & S. Gupta, EIR-Bericht 247, November 1973
C              Eidg. Institut fur Reaktorforschung
C              Wuerenlingen, Switzerland
C   see also:  Haenggi, Roesel & Trautmann,
C              J. Comput. Phys., v. 137, (1980) 242-258
C   note:      restart procedure modified by I.J.Thompson
C
C*******************************************************************
C
      IMPLICIT COMPLEX*16(A-H,O-Z)
      DIMENSION A(100),B(100),XX(2,100)
      LOGICAL EVEN
      DOUBLE PRECISION EPS

      IBN=INUM
      IF(IBEG .GT. 4) GO TO 50
      IF(IBEG .EQ. 4) GO TO 20
      B(1)=A(1)
      IF(IBN .GE. 2) B(2)=-A(2)/A(1)
      IF(IBN .LT. 3) RETURN
      X0=A(3)/A(2)
      XX(2,1)=B(2)
      XX(1,1)=-X0
      XX(1,2)=0
      B(3)=-X0-B(2)
      X0=-B(3)*A(2)
      M=3
      MP12=2
      EVEN=.TRUE.
      IF(IBN .LE. 3) RETURN
   20 IF(ABS(B(3)) .LT. EPS*ABS(X0)) THEN
       INUM=-M
       RETURN
      END IF
      M=4
   30 X1=A(M)
      M2M1=MP12
      MP12=M2M1+1
      IF(EVEN) MP12=M2M1
      DO 40 K = 2,MP12
   40 X1=X1+A(M-K+1)*XX(1,K-1)
      B(M)=-X1/X0
      IF(M .GE. IBN) RETURN
   50 IF(ABS(B(M)) .LT. EPS*ABS(X0)) THEN
       INUM=-M
       RETURN
      END IF
      DO 60 K = M2M1,2,-1
   60 XX(2,K)=XX(1,K)+B(M)*XX(2,K-1)
      XX(2,1)=XX(1,1)+B(M)
      DO 70 K = 1,M2M1
      X0=XX(2,K)
      XX(2,K)=XX(1,K)
   70 XX(1,K)=X0
      X0=X1
      XX(1,M2M1+1)=0
      M=M+1
      EVEN=.NOT.EVEN
      GO TO 30
      END
#endif
Commit	Line	Data
fe4da5cc	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/z3 + A3/z5 + ... + 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