C	as dpc.for but for a linear system and with scaling
C	modification by Sven Ivansson  1990



	SUBROUTINE DPCLIN ( FSUBLIN,N, EPSGLOBLN,EPSGLOB,EPSFAKT,
     $	HFAKTMAX,HOFF,EPSMAGNIFMAX, XEND, X0,DPCLN,Y0,EPS,HSTEP, IND )
C	Dormand-Prince Complex for a LINear system
C	scaling is used
C	this routine solves a system of first-order ordinary
C	differential equations by Dormand-Prince of local order 5 (4) 
C	six function evaluations per step
C	taken from  Hairer,Noersett,Wanner (1987) "Solving Ordinary 
C	Differential Equations: I. Non-stiff problems", p.434
C	    y' = f(x,y) ,  y(x0) = y0
C	input parameters:
C	   fsublin         external subroutine :    fsublin ( n,x,y, dy )
C	                      input:  integer n
C	                              real*8 x
C	                              complex*16 y(1:n)
C	                      output: complex*16 dy(1:n)   " dy = y'(x) "
C	                      NOTE:  y' is assumed to be linear in y
C	   n               number of equations
C	   epsglobln       global error tolerance for dpcln  (absolute error, 
C	                   but actually more like relative error because of 
C	                   the scaling)
C	                   (when epsglobln is nonpositive, it is assumed to be 
C	                   infinite)
C	   epsglob         global error tolerance for the scaled y0 vector 
C	                   (absolute error, but actually more like relative 
C	                   error because of the scaling)
C	                   (when epsglob is nonpositive, it is assumed to be 
C	                   infinite)
C	   epsfakt         reduction factor for eps  (used only when epsglobln 
C	                   or epsglob is positive)
C	                   ( take  0.0 < epsfakt < 1.0  )
C	   hfaktmax        maximum magnification factor for dabs(hlok) in 
C	                   each step          ( 1.05d0 might be suitable )
C	   hoff            lower limit for dabs(hlok)   ( take  hoff >= 0.0d0 )
C	   epsmagnifmax    maximum value for magnified output eps
C	   xend            final x-value ("both directions" possible)
C	input/output parameters:
C	   x0              initial/final x-value
C	   dpcln           initial/final exponential de-scaling for y0
C	   y0              initial/final y-value is   dexp(dpcln) * y0
C	                   initially, it is automatically arranged that 
C	                     dmax1(dabs(y0(1:n)) := "1.0d0 or 0.0d0"  with the 
C	                   implied change of dpcln
C	   eps             "first turn"/"last (but zero, one, or two) turn" 
C	                   local error tolerance for the "one-sided scaled" 
C	                   y0-vector     (absolute error, but actually more 
C	                   like relative error because of the scaling)
C	                   (when eps is nonpositive, it is assumed to be 
C	                   infinite; positive epsglobln or epsglob are 
C	                   then not allowed, hfaktmax and hoff and 
C	                   epsmagnifmax are then irrelevant except that 
C	                   hoff must not be greater than hstep)
C	                   ( eps = 1d-8 on entry might be suitable ) 
C	                   ( eps will be unchanged if  epsglob <= 0.0d0 & 
C	                   epsglobln <= 0.0d0 )
C	   hstep           "first turn initial"/"last (but zero or one) turn 
C	                   final-proposal-upon-acceptance" absolute step size 
C	                   (for a periodic function, it might be suitable to 
C	                   take "hstep on entry" as one quarter of a period)
C	                   ( hstep will be unchanged if  eps <= 0.0d0 )
C	output parameters:
C	   ind = 0         successful return with  x0 = xend
C	       = 1         unsuccessful return   ( dabs(hlok) too small )
C	output "common":
C	   /dpcstat/

	IMPLICIT NONE

	INTEGER NDIM
	REAL*8 FAC1,FAC2, ERREXP, CONTRUNC, HFAKTZERO, 
     $	DMAXLIM1,DMAXLIM2
	PARAMETER ( NDIM = 100 , FAC1 = 0.9D0 , FAC2 = 0.9D0 , 
     $	ERREXP = 0.2D0 , CONTRUNC = 1.0D-10 , HFAKTZERO = 1.0D-5 , 
     $	DMAXLIM1 = 0.01D0 , DMAXLIM2 = 100.0D0 )

	INTEGER N
	REAL*8 EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,EPSMAGNIFMAX, 
     $	XEND
	EXTERNAL FSUBLIN

	REAL*8 X0, DPCLN, EPS,HSTEP
	COMPLEX*16 Y0(*)

	INTEGER IND

	COMMON /DPCSTAT/ NDPCTURN,NDPCACCPT,NDPCREJCT
	INTEGER NDPCTURN,NDPCACCPT,NDPCREJCT

	INTEGER I
	REAL*8 RCA1,RCA2,RCA3, RCB1,RCB2,RCB3,RCB4, RCC1,RCC2,RCC3,
     $	RCC4,RCC5, RCD1,RCD2,RCD3,RCD4,RCD5, RCE1,RCE2,RCE3,RCE4,RCE5, 
     $	SRIKT, X0INIT, X,X1, DPCLNINIT, DPCLNOLD, EPSOLD,HSTEPOLD, TL, 
     $	HLOK, HABS,H5, CON,ERR, SLASKLN, SLASK
	COMPLEX*16 Y0INIT(ndim), Y0OLD(ndim), Y1(ndim), 
     $	CK1(ndim), CK2(ndim), CK3(ndim), CK4(ndim), CK5(ndim)
	LOGICAL FIRST

	SAVE RCA1,RCA2,RCA3, RCB1,RCB2,RCB3,RCB4, RCC1,RCC2,RCC3,
     $	RCC4,RCC5, RCD1,RCD2,RCD3,RCD4,RCD5, RCE1,RCE2,RCE3,RCE4,RCE5
	DATA FIRST /.TRUE./

	IF (FIRST) THEN

	  RCA1 = 44.0D0 / 45.0D0
	  RCA2 = 56.0D0 / 15.0D0
	  RCA3 = 32.0D0 / 9.0D0

	  RCB1 = 19372.0D0 / 6561.0D0
	  RCB2 = 25360.0D0 / 2187.0D0
	  RCB3 = 64448.0D0 / 6561.0D0
	  RCB4 = 212.0D0 / 729.0D0

	  RCC1 = 9017.0D0 / 3168.0D0
	  RCC2 = 355.0D0 / 33.0D0
	  RCC3 = 46732.0D0 / 5247.0D0
	  RCC4 = 49.0D0 / 176.0D0
	  RCC5 = 5103.0D0 / 18656.0D0

	  RCD1 = 35.0D0 / 384.0D0
	  RCD2 = 500.0D0 / 1113.0D0
	  RCD3 = 125.0D0 / 192.0D0
	  RCD4 = 2187.0D0 / 6784.0D0
	  RCD5 = 11.0D0 / 84.0D0

	  RCE1 = 71.0D0 / 57600.0D0
	  RCE2 = 71.0D0 / 16695.0D0
	  RCE3 = 71.0D0 / 1920.0D0
	  RCE4 = 17253.0D0 / 339200.0D0
	  RCE5 = 22.0D0 / 525.0D0

	  FIRST = .FALSE.

	END IF

	IF (N.GT.ndim) STOP 540

	NDPCTURN = 0
	NDPCACCPT = 0
	NDPCREJCT = 0

	IF (X0.EQ.XEND) RETURN
	IF (HSTEP.LE.HOFF) THEN
	  IND = 1
	  RETURN
	END IF

	SLASK = 0.0D0
	DO I = 1, N
	  SLASK = DMAX1 ( SLASK , DABS(DREAL(Y0(I))), DABS(DIMAG(Y0(I))) )
	END DO
	IF (SLASK.EQ.0.0D0) RETURN
	IF ( SLASK.LT.dmaxlim1 .OR. dmaxlim2.LT.SLASK ) THEN
	  SLASK = 1/SLASK
	  DPCLN = DPCLN - DLOG(SLASK)
	  DO I = 1, N
	    Y0(I) = SLASK * Y0(I)
	  END DO
	END IF

	SRIKT = DSIGN ( 1.0D0 , XEND-X0 )

	IF ( EPSGLOBLN.GT.0.0D0 .OR.  EPSGLOB.GT.0.0D0 ) THEN
	  IF ( EPSFAKT.LE.0.0D0 .OR. 1.0D0.LE.EPSFAKT ) STOP 460
	  IF (EPS.LE.0.0D0) STOP 461
	  X0INIT = X0
	  DPCLNINIT = DPCLN
	  DO I = 1, N
	    Y0INIT(I) = Y0(I)
	  END DO
	END IF

	DO WHILE (.TRUE.)

	  NDPCTURN = NDPCTURN + 1
	  IF (NDPCTURN.GT.1) THEN
	    X0 = X0INIT
	    DPCLNOLD = DPCLN
	    DPCLN = DPCLNINIT
	    DO I = 1, N
	      Y0OLD(I) = Y0(I)
	      Y0(I) = Y0INIT(I)
	    END DO
	    EPSOLD = EPS
	    HSTEPOLD = HSTEP
	    EPS = EPS * EPSFAKT
	  END IF

	  IF (EPS.GT.0.0D0) TL = EPS**ERREXP

	  CALL FSUBLIN ( N,X0,Y0, CK1 )

	  HLOK = DSIGN (HSTEP,SRIKT)

	  DO WHILE (X0.NE.XEND)

	    X1 = X0 + HLOK
	    IF (X1.EQ.X0) THEN
	      IND = 1
	      RETURN
	    END IF
	    IF (SRIKT.GT.0.0D0) THEN
	      IF (X1.GT.XEND) THEN
	        X1 = XEND
	        HLOK = XEND - X0
	      END IF
	    ELSE
	      IF (X1.LT.XEND) THEN
	        X1 = XEND
	        HLOK = XEND - X0
	      END IF
	    END IF
	    HABS = DABS(HLOK)
	    H5 = HABS**5

	    SLASK = 0.2D0 * HLOK
	    DO I = 1, N
	      Y1(I) = Y0(I) + SLASK * CK1(I)
	    END DO

	    X = X0 + 0.2D0*HLOK
	    CALL FSUBLIN ( N,X,Y1, CK2 )
	    DO I = 1, N
	      Y1(I) = Y0(I) + HLOK * ( 0.075D0*CK1(I) + 0.225D0*CK2(I) )
	    END DO

	    X = X0 + 0.3D0*HLOK
	    CALL FSUBLIN ( N,X,Y1, CK3 )
	    DO I = 1, N
	      Y1(I) = Y0(I) + HLOK * ( RCA1*CK1(I) - RCA2*CK2(I) + 
     $	      RCA3*CK3(I) )
	    END DO

	    X = X0 + 0.8D0*HLOK
	    CALL FSUBLIN ( N,X,Y1, CK4 )
	    DO I = 1, N
	      Y1(I) = Y0(I) + HLOK * ( RCB1*CK1(I) - RCB2*CK2(I) + 
     $	      RCB3*CK3(I) - RCB4*CK4(I) )
	    END DO

	    X = X0 + 8.0D0*HLOK/9.0D0
	    CALL FSUBLIN ( N,X,Y1, CK5 )
	    DO I = 1, N
	      Y1(I) = Y0(I) + HLOK * ( RCC1*CK1(I) - RCC2*CK2(I) + 
     $	      RCC3*CK3(I) + RCC4*CK4(I) - RCC5*CK5(I) )
	    END DO

	    CALL FSUBLIN ( N,X1,Y1, CK2 )
	    DO I = 1, N
	      Y1(I) = Y0(I) + HLOK * ( RCD1*CK1(I) + RCD2*CK3(I) + 
     $	      RCD3*CK4(I) - RCD4*CK5(I) + RCD5*CK2(I))
	    END DO

C	    compute intermediate sum to save memory
	    DO I = 1, N
	      CK2(I) = RCE1*CK1(I) - RCE2*CK3(I) + RCE3*CK4(I) - 
     $	      RCE4*CK5(I) + RCE5*CK2(I)
	    END DO

C 	    the last stage
	    CALL FSUBLIN ( N,X1,Y1, CK3 )
	    ERR = 0.0D0
	    DO I = 1, N
	      CK4(I) = ( CK2(I)-0.025D0*CK3(I) ) * HLOK
	      ERR = DMAX1 ( ERR , DABS(DREAL(CK4(I))), DABS(DIMAG(CK4(I))) )
	    END DO

	    IF ( EPS.LE.0.0D0 .OR. ERR.LE.EPS ) THEN
C	      step is accepted
	      NDPCACCPT = NDPCACCPT + 1
	      IF (EPS.GT.0.0D0) THEN
	        IF (H5.GT.0.0D0) THEN
	          CON = DMAX1 ( CONTRUNC , ERR/H5 )
	          CON  = CON**ERREXP
	          HSTEP = FAC1*TL/CON
	        ELSE
	          HSTEP = HFAKTZERO*HABS
	        END IF
	        HLOK = DMIN1 ( HSTEP , HFAKTMAX*HABS )
	        IF (HLOK.LE.HOFF) THEN
	          IND = 1
	          RETURN
	        END IF
	        HLOK = DSIGN (HLOK,SRIKT)
	      END IF
C	      (do the scaling now)
	      SLASK = 0.0D0
	      DO I = 1, N
	        SLASK = DMAX1 ( SLASK , DABS(DREAL(Y1(I))), 
     $	        DABS(DIMAG(Y1(I))) )
	      END DO
	      IF ( SLASK.LT.dmaxlim1 .OR. dmaxlim2.LT.SLASK ) THEN
	        SLASK = 1/SLASK
	        DPCLN = DPCLN - DLOG(SLASK)
	        DO I = 1, N
	          CK1(I) = SLASK * CK3(I)
	          Y0(I) = SLASK * Y1(I)
	        END DO
	      ELSE
	        DO I = 1, N
	          CK1(I) = CK3(I)
	          Y0(I) = Y1(I)
	        END DO
	      END IF
	      X0 = X1
	    ELSE
C	      step is rejected
	      NDPCREJCT = NDPCREJCT + 1
	      HLOK = FAC2 * TL * HABS / ERR**ERREXP
	      IF (HLOK.LE.HOFF) THEN
	        IND = 1
	        RETURN
	      END IF
	      HLOK = DSIGN (HLOK,SRIKT)
	    END IF

	  END DO

	  IF (NDPCTURN.EQ.1) THEN
	    IF ( EPSGLOBLN.LE.0.0D0 .AND. EPSGLOB.LE.0.0D0 ) THEN
	      IND = 0
	      RETURN
	    END IF
	  ELSE
	    IF (EPSGLOBLN.GT.0.0D0) THEN
	      SLASKLN = DABS(DPCLN-DPCLNOLD)
	    ELSE
	      SLASKLN = EPSGLOBLN
	    END IF
	    IF (EPSGLOB.GT.0.0D0) THEN
	      SLASK = 0.0D0
	      DO I = 1, N
	        SLASK = DMAX1 ( SLASK , DABS(DREAL(Y0OLD(I)-Y0(I))), 
     $	        DABS(DIMAG(Y0OLD(I)-Y0(I))) )
	      END DO
	    ELSE
	      SLASK = EPSGLOB
	    END IF
	    IF ( SLASKLN.LE.EPSGLOBLN .AND. SLASK.LE.EPSGLOB ) THEN
	      EPS = EPSOLD
	      HSTEP = HSTEPOLD
	      IF (NDPCTURN.EQ.2) THEN
	        IF (SLASKLN.LE.(EPSFAKT*EPSGLOBLN)) THEN
	          IF (SLASK.LE.(EPSFAKT*EPSGLOB)) THEN
	            EPS =  DMIN1 ( EPS/EPSFAKT , EPSMAGNIFMAX )
	          END IF
	        END IF
	      END IF
	      IND = 0
	      RETURN
	    END IF
	  END IF

	END DO

	END