C	Ilkka Karasalo's EULETRAC.FOR mildly adapted by Sven Ivansson 1990
 
C---------------------------------------------------------------
C
C	EULERACC
C	computes the sum of a COMPLEX*16 valued series with slowly 
C	decreasing, oscillating terms 
C	convergence acceleration by the Euler transformation
C
C	Euler's transformation (see DBA, Sec 7.6, Problem 3), is intended 
C	for a series whose n:th term is of the form 
C	          term(n)  =  u(n) * (zosc**n)
C	where  cdabs(zosc) <= 1  (preferably zosc should be close to -1)
C	and  u(n) is smooth
C
C---------------------------------------------------------------

	SUBROUTINE EULERACC (MTOT,TERFUN,ZOSC, NSTA,NEND, KMAX,EPS, 
     $	N,K,SUM, ICWORKLIM0,ICWORKDIM,CWORK)
C	mtot         dimension of vector-valued terms
C	terfun       name of routine such that  
C	             " call terfun (n, term, icworklim0,icworkdim,cwork) "  
C	             returns for given n with  nsta <= n <= nend  
C	                term(1:mtot)    n:th (vector-valued) term in series
C	zosc         base of (the oscillating) geometric factors of the terms
C	nsta,nend    starting value and final value (unless convergence 
C	             occurs before final value is reached) of index n of 
C	             the series,  "nend < nsta"  implies  "nend = +infinity"
C	kmax         max number of columns to be used in the difference 
C	             scheme = 1 + max number of terms in the transformed 
C	             series,  must have  2 <= kmax
C	eps          absolute tolerance (strictly valid under certain 
C	             conditions, see DBA)
C	n            index of last term needed (value of n at last call to 
C	             "terfun")  (output)
C	k            column index in difference table for the accepted 
C	             elements  (output)
C	sum          (mtot-dimensional vector-valued) sum  (output)
C	icworklim0
C	icworkdim
C	cwork(1:icworkdim)   "containing cwork(icworklim0+1,icworkdim)" 
C	                     for user's own work-space  < but note that 
C	                     cwork(1:icworklim0) must not be accessed here >

	IMPLICIT NONE

	INTEGER MTOTDIM, KMAXDIM
	PARAMETER ( MTOTDIM = 550 , KMAXDIM = 40 )

	INTEGER NSTA,NEND, KMAX
	REAL*8 EPS
	COMPLEX*16 ZOSC
	EXTERNAL TERFUN

	INTEGER N, K
	COMPLEX*16 SUM(*)

	INTEGER ICWORKLIM0,ICWORKDIM
	COMPLEX*16 CWORK(*)

	INTEGER MTOT, KCOL, I, KM1
	REAL*8 ZOSCA
	COMPLEX*16 ZEXP, CSLASK, 
     $	TERM(mtotdim), TABLE(mtotdim,kmaxdim), DIFF(mtotdim), 
     $	VOLD(mtotdim)

	IF (MTOT.GT.mtotdim .OR. KMAX.LT.2 .OR. 
     $	KMAX.GT.(kmaxdim+1) ) STOP 111
 
	ZOSCA = CDABS(ZOSC)
	ZEXP = 1.0D0 - ZOSC
	IF (ZEXP.EQ.DCMPLX(0.D0,0.D0)) STOP 222
	ZEXP = 1.0D0/ZEXP     ! zexp := 1/(1-zosc)

	N = NSTA
	CALL TERFUN (N, TERM, ICWORKLIM0,ICWORKDIM,CWORK)
	DO I = 1, MTOT
	  TABLE(I,1) = TERM(I)
	END DO

	KCOL = 2
 
	DO WHILE (N.NE.NEND)    ! loop over terms

	  N = N + 1
	  CALL TERFUN (N, TERM, ICWORKLIM0,ICWORKDIM,CWORK)
	  DO I = 1, MTOT
	    VOLD(I) = TABLE(I,1)
	    TABLE(I,1) = VOLD(I) + TERM(I)
	  END DO
 
	  DO K = 2, KCOL    ! loop over columns
 
	    KM1 = K - 1
	    DO I = 1, MTOT
	      DIFF(I) = ZEXP * (TABLE(I,KM1)-VOLD(I))  ! * "vert diff in col k"
	    END DO
 
	    DO I = 1, MTOT
	      IF ( DMAX1(DABS(DREAL(DIFF(I))),DABS(DIMAG(DIFF(I))))*
     $	      ZOSCA .GT. EPS ) GOTO 11  ! convergence test
	    END DO
	    DO I = 1, MTOT
	      SUM(I) = VOLD(I) + DIFF(I)
	    END DO
	    RETURN

 11	    IF (K.LT.KMAX) THEN
	      DO I = 1, MTOT
	        CSLASK = VOLD(I)
	        VOLD(I) = TABLE(I,K)
	        TABLE(I,K) = CSLASK + DIFF(I)     ! update next column
	      END DO
	    END IF
 
	  END DO   ! end extrapolation loop
 
	  IF (KCOL.LT.KMAX) KCOL = KCOL + 1

	END DO  ! end loop over terms
 
	STOP 66      ! no convergence

	END