[Scipy-svn] r3951 - in trunk/scipy/integrate: linpack_lite odepack
scipy-svn@scip...
scipy-svn@scip...
Tue Feb 19 23:39:10 CST 2008
Author: oliphant
Date: 2008-02-19 23:39:03 -0600 (Tue, 19 Feb 2008)
New Revision: 3951
Added:
trunk/scipy/integrate/linpack_lite/zgbfa.f
trunk/scipy/integrate/linpack_lite/zgbsl.f
trunk/scipy/integrate/linpack_lite/zgefa.f
trunk/scipy/integrate/linpack_lite/zgesl.f
trunk/scipy/integrate/odepack/zvode.f
Log:
Add missing pieces from ticket #334
Added: trunk/scipy/integrate/linpack_lite/zgbfa.f
===================================================================
--- trunk/scipy/integrate/linpack_lite/zgbfa.f 2008-02-20 05:27:16 UTC (rev 3950)
+++ trunk/scipy/integrate/linpack_lite/zgbfa.f 2008-02-20 05:39:03 UTC (rev 3951)
@@ -0,0 +1,181 @@
+ subroutine zgbfa(abd,lda,n,ml,mu,ipvt,info)
+ integer lda,n,ml,mu,ipvt(1),info
+ complex*16 abd(lda,1)
+c
+c zgbfa factors a complex*16 band matrix by elimination.
+c
+c zgbfa is usually called by zgbco, but it can be called
+c directly with a saving in time if rcond is not needed.
+c
+c on entry
+c
+c abd complex*16(lda, n)
+c contains the matrix in band storage. the columns
+c of the matrix are stored in the columns of abd and
+c the diagonals of the matrix are stored in rows
+c ml+1 through 2*ml+mu+1 of abd .
+c see the comments below for details.
+c
+c lda integer
+c the leading dimension of the array abd .
+c lda must be .ge. 2*ml + mu + 1 .
+c
+c n integer
+c the order of the original matrix.
+c
+c ml integer
+c number of diagonals below the main diagonal.
+c 0 .le. ml .lt. n .
+c
+c mu integer
+c number of diagonals above the main diagonal.
+c 0 .le. mu .lt. n .
+c more efficient if ml .le. mu .
+c on return
+c
+c abd an upper triangular matrix in band storage and
+c the multipliers which were used to obtain it.
+c the factorization can be written a = l*u where
+c l is a product of permutation and unit lower
+c triangular matrices and u is upper triangular.
+c
+c ipvt integer(n)
+c an integer vector of pivot indices.
+c
+c info integer
+c = 0 normal value.
+c = k if u(k,k) .eq. 0.0 . this is not an error
+c condition for this subroutine, but it does
+c indicate that zgbsl will divide by zero if
+c called. use rcond in zgbco for a reliable
+c indication of singularity.
+c
+c band storage
+c
+c if a is a band matrix, the following program segment
+c will set up the input.
+c
+c ml = (band width below the diagonal)
+c mu = (band width above the diagonal)
+c m = ml + mu + 1
+c do 20 j = 1, n
+c i1 = max0(1, j-mu)
+c i2 = min0(n, j+ml)
+c do 10 i = i1, i2
+c k = i - j + m
+c abd(k,j) = a(i,j)
+c 10 continue
+c 20 continue
+c
+c this uses rows ml+1 through 2*ml+mu+1 of abd .
+c in addition, the first ml rows in abd are used for
+c elements generated during the triangularization.
+c the total number of rows needed in abd is 2*ml+mu+1 .
+c the ml+mu by ml+mu upper left triangle and the
+c ml by ml lower right triangle are not referenced.
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas zaxpy,zscal,izamax
+c fortran dabs,max0,min0
+c
+c internal variables
+c
+ complex*16 t
+ integer i,izamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
+c
+ complex*16 zdum
+ double precision cabs1
+ double precision dreal,dimag
+ complex*16 zdumr,zdumi
+ dreal(zdumr) = zdumr
+ dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
+ cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
+c
+ m = ml + mu + 1
+ info = 0
+c
+c zero initial fill-in columns
+c
+ j0 = mu + 2
+ j1 = min0(n,m) - 1
+ if (j1 .lt. j0) go to 30
+ do 20 jz = j0, j1
+ i0 = m + 1 - jz
+ do 10 i = i0, ml
+ abd(i,jz) = (0.0d0,0.0d0)
+ 10 continue
+ 20 continue
+ 30 continue
+ jz = j1
+ ju = 0
+c
+c gaussian elimination with partial pivoting
+c
+ nm1 = n - 1
+ if (nm1 .lt. 1) go to 130
+ do 120 k = 1, nm1
+ kp1 = k + 1
+c
+c zero next fill-in column
+c
+ jz = jz + 1
+ if (jz .gt. n) go to 50
+ if (ml .lt. 1) go to 50
+ do 40 i = 1, ml
+ abd(i,jz) = (0.0d0,0.0d0)
+ 40 continue
+ 50 continue
+c
+c find l = pivot index
+c
+ lm = min0(ml,n-k)
+ l = izamax(lm+1,abd(m,k),1) + m - 1
+ ipvt(k) = l + k - m
+c
+c zero pivot implies this column already triangularized
+c
+ if (cabs1(abd(l,k)) .eq. 0.0d0) go to 100
+c
+c interchange if necessary
+c
+ if (l .eq. m) go to 60
+ t = abd(l,k)
+ abd(l,k) = abd(m,k)
+ abd(m,k) = t
+ 60 continue
+c
+c compute multipliers
+c
+ t = -(1.0d0,0.0d0)/abd(m,k)
+ call zscal(lm,t,abd(m+1,k),1)
+c
+c row elimination with column indexing
+c
+ ju = min0(max0(ju,mu+ipvt(k)),n)
+ mm = m
+ if (ju .lt. kp1) go to 90
+ do 80 j = kp1, ju
+ l = l - 1
+ mm = mm - 1
+ t = abd(l,j)
+ if (l .eq. mm) go to 70
+ abd(l,j) = abd(mm,j)
+ abd(mm,j) = t
+ 70 continue
+ call zaxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
+ 80 continue
+ 90 continue
+ go to 110
+ 100 continue
+ info = k
+ 110 continue
+ 120 continue
+ 130 continue
+ ipvt(n) = n
+ if (cabs1(abd(m,n)) .eq. 0.0d0) info = n
+ return
+ end
Property changes on: trunk/scipy/integrate/linpack_lite/zgbfa.f
___________________________________________________________________
Name: svn:eol-style
+ native
Added: trunk/scipy/integrate/linpack_lite/zgbsl.f
===================================================================
--- trunk/scipy/integrate/linpack_lite/zgbsl.f 2008-02-20 05:27:16 UTC (rev 3950)
+++ trunk/scipy/integrate/linpack_lite/zgbsl.f 2008-02-20 05:39:03 UTC (rev 3951)
@@ -0,0 +1,139 @@
+ subroutine zgbsl(abd,lda,n,ml,mu,ipvt,b,job)
+ integer lda,n,ml,mu,ipvt(1),job
+ complex*16 abd(lda,1),b(1)
+c
+c zgbsl solves the complex*16 band system
+c a * x = b or ctrans(a) * x = b
+c using the factors computed by zgbco or zgbfa.
+c
+c on entry
+c
+c abd complex*16(lda, n)
+c the output from zgbco or zgbfa.
+c
+c lda integer
+c the leading dimension of the array abd .
+c
+c n integer
+c the order of the original matrix.
+c
+c ml integer
+c number of diagonals below the main diagonal.
+c
+c mu integer
+c number of diagonals above the main diagonal.
+c
+c ipvt integer(n)
+c the pivot vector from zgbco or zgbfa.
+c
+c b complex*16(n)
+c the right hand side vector.
+c
+c job integer
+c = 0 to solve a*x = b ,
+c = nonzero to solve ctrans(a)*x = b , where
+c ctrans(a) is the conjugate transpose.
+c
+c on return
+c
+c b the solution vector x .
+c
+c error condition
+c
+c a division by zero will occur if the input factor contains a
+c zero on the diagonal. technically this indicates singularity
+c but it is often caused by improper arguments or improper
+c setting of lda . it will not occur if the subroutines are
+c called correctly and if zgbco has set rcond .gt. 0.0
+c or zgbfa has set info .eq. 0 .
+c
+c to compute inverse(a) * c where c is a matrix
+c with p columns
+c call zgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
+c if (rcond is too small) go to ...
+c do 10 j = 1, p
+c call zgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
+c 10 continue
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas zaxpy,zdotc
+c fortran dconjg,min0
+c
+c internal variables
+c
+ complex*16 zdotc,t
+ integer k,kb,l,la,lb,lm,m,nm1
+ double precision dreal,dimag
+ complex*16 zdumr,zdumi
+ dreal(zdumr) = zdumr
+ dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
+c
+ m = mu + ml + 1
+ nm1 = n - 1
+ if (job .ne. 0) go to 50
+c
+c job = 0 , solve a * x = b
+c first solve l*y = b
+c
+ if (ml .eq. 0) go to 30
+ if (nm1 .lt. 1) go to 30
+ do 20 k = 1, nm1
+ lm = min0(ml,n-k)
+ l = ipvt(k)
+ t = b(l)
+ if (l .eq. k) go to 10
+ b(l) = b(k)
+ b(k) = t
+ 10 continue
+ call zaxpy(lm,t,abd(m+1,k),1,b(k+1),1)
+ 20 continue
+ 30 continue
+c
+c now solve u*x = y
+c
+ do 40 kb = 1, n
+ k = n + 1 - kb
+ b(k) = b(k)/abd(m,k)
+ lm = min0(k,m) - 1
+ la = m - lm
+ lb = k - lm
+ t = -b(k)
+ call zaxpy(lm,t,abd(la,k),1,b(lb),1)
+ 40 continue
+ go to 100
+ 50 continue
+c
+c job = nonzero, solve ctrans(a) * x = b
+c first solve ctrans(u)*y = b
+c
+ do 60 k = 1, n
+ lm = min0(k,m) - 1
+ la = m - lm
+ lb = k - lm
+ t = zdotc(lm,abd(la,k),1,b(lb),1)
+ b(k) = (b(k) - t)/dconjg(abd(m,k))
+ 60 continue
+c
+c now solve ctrans(l)*x = y
+c
+ if (ml .eq. 0) go to 90
+ if (nm1 .lt. 1) go to 90
+ do 80 kb = 1, nm1
+ k = n - kb
+ lm = min0(ml,n-k)
+ b(k) = b(k) + zdotc(lm,abd(m+1,k),1,b(k+1),1)
+ l = ipvt(k)
+ if (l .eq. k) go to 70
+ t = b(l)
+ b(l) = b(k)
+ b(k) = t
+ 70 continue
+ 80 continue
+ 90 continue
+ 100 continue
+ return
+ end
Property changes on: trunk/scipy/integrate/linpack_lite/zgbsl.f
___________________________________________________________________
Name: svn:eol-style
+ native
Added: trunk/scipy/integrate/linpack_lite/zgefa.f
===================================================================
--- trunk/scipy/integrate/linpack_lite/zgefa.f 2008-02-20 05:27:16 UTC (rev 3950)
+++ trunk/scipy/integrate/linpack_lite/zgefa.f 2008-02-20 05:39:03 UTC (rev 3951)
@@ -0,0 +1,111 @@
+ subroutine zgefa(a,lda,n,ipvt,info)
+ integer lda,n,ipvt(1),info
+ complex*16 a(lda,1)
+c
+c zgefa factors a complex*16 matrix by gaussian elimination.
+c
+c zgefa is usually called by zgeco, but it can be called
+c directly with a saving in time if rcond is not needed.
+c (time for zgeco) = (1 + 9/n)*(time for zgefa) .
+c
+c on entry
+c
+c a complex*16(lda, n)
+c the matrix to be factored.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c on return
+c
+c a an upper triangular matrix and the multipliers
+c which were used to obtain it.
+c the factorization can be written a = l*u where
+c l is a product of permutation and unit lower
+c triangular matrices and u is upper triangular.
+c
+c ipvt integer(n)
+c an integer vector of pivot indices.
+c
+c info integer
+c = 0 normal value.
+c = k if u(k,k) .eq. 0.0 . this is not an error
+c condition for this subroutine, but it does
+c indicate that zgesl or zgedi will divide by zero
+c if called. use rcond in zgeco for a reliable
+c indication of singularity.
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas zaxpy,zscal,izamax
+c fortran dabs
+c
+c internal variables
+c
+ complex*16 t
+ integer izamax,j,k,kp1,l,nm1
+c
+ complex*16 zdum
+ double precision cabs1
+ double precision dreal,dimag
+ complex*16 zdumr,zdumi
+ dreal(zdumr) = zdumr
+ dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
+ cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
+c
+c gaussian elimination with partial pivoting
+c
+ info = 0
+ nm1 = n - 1
+ if (nm1 .lt. 1) go to 70
+ do 60 k = 1, nm1
+ kp1 = k + 1
+c
+c find l = pivot index
+c
+ l = izamax(n-k+1,a(k,k),1) + k - 1
+ ipvt(k) = l
+c
+c zero pivot implies this column already triangularized
+c
+ if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
+c
+c interchange if necessary
+c
+ if (l .eq. k) go to 10
+ t = a(l,k)
+ a(l,k) = a(k,k)
+ a(k,k) = t
+ 10 continue
+c
+c compute multipliers
+c
+ t = -(1.0d0,0.0d0)/a(k,k)
+ call zscal(n-k,t,a(k+1,k),1)
+c
+c row elimination with column indexing
+c
+ do 30 j = kp1, n
+ t = a(l,j)
+ if (l .eq. k) go to 20
+ a(l,j) = a(k,j)
+ a(k,j) = t
+ 20 continue
+ call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
+ 30 continue
+ go to 50
+ 40 continue
+ info = k
+ 50 continue
+ 60 continue
+ 70 continue
+ ipvt(n) = n
+ if (cabs1(a(n,n)) .eq. 0.0d0) info = n
+ return
+ end
Property changes on: trunk/scipy/integrate/linpack_lite/zgefa.f
___________________________________________________________________
Name: svn:eol-style
+ native
Added: trunk/scipy/integrate/linpack_lite/zgesl.f
===================================================================
--- trunk/scipy/integrate/linpack_lite/zgesl.f 2008-02-20 05:27:16 UTC (rev 3950)
+++ trunk/scipy/integrate/linpack_lite/zgesl.f 2008-02-20 05:39:03 UTC (rev 3951)
@@ -0,0 +1,122 @@
+ subroutine zgesl(a,lda,n,ipvt,b,job)
+ integer lda,n,ipvt(1),job
+ complex*16 a(lda,1),b(1)
+c
+c zgesl solves the complex*16 system
+c a * x = b or ctrans(a) * x = b
+c using the factors computed by zgeco or zgefa.
+c
+c on entry
+c
+c a complex*16(lda, n)
+c the output from zgeco or zgefa.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c ipvt integer(n)
+c the pivot vector from zgeco or zgefa.
+c
+c b complex*16(n)
+c the right hand side vector.
+c
+c job integer
+c = 0 to solve a*x = b ,
+c = nonzero to solve ctrans(a)*x = b where
+c ctrans(a) is the conjugate transpose.
+c
+c on return
+c
+c b the solution vector x .
+c
+c error condition
+c
+c a division by zero will occur if the input factor contains a
+c zero on the diagonal. technically this indicates singularity
+c but it is often caused by improper arguments or improper
+c setting of lda . it will not occur if the subroutines are
+c called correctly and if zgeco has set rcond .gt. 0.0
+c or zgefa has set info .eq. 0 .
+c
+c to compute inverse(a) * c where c is a matrix
+c with p columns
+c call zgeco(a,lda,n,ipvt,rcond,z)
+c if (rcond is too small) go to ...
+c do 10 j = 1, p
+c call zgesl(a,lda,n,ipvt,c(1,j),0)
+c 10 continue
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas zaxpy,zdotc
+c fortran dconjg
+c
+c internal variables
+c
+ complex*16 zdotc,t
+ integer k,kb,l,nm1
+ double precision dreal,dimag
+ complex*16 zdumr,zdumi
+ dreal(zdumr) = zdumr
+ dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
+c
+ nm1 = n - 1
+ if (job .ne. 0) go to 50
+c
+c job = 0 , solve a * x = b
+c first solve l*y = b
+c
+ if (nm1 .lt. 1) go to 30
+ do 20 k = 1, nm1
+ l = ipvt(k)
+ t = b(l)
+ if (l .eq. k) go to 10
+ b(l) = b(k)
+ b(k) = t
+ 10 continue
+ call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1)
+ 20 continue
+ 30 continue
+c
+c now solve u*x = y
+c
+ do 40 kb = 1, n
+ k = n + 1 - kb
+ b(k) = b(k)/a(k,k)
+ t = -b(k)
+ call zaxpy(k-1,t,a(1,k),1,b(1),1)
+ 40 continue
+ go to 100
+ 50 continue
+c
+c job = nonzero, solve ctrans(a) * x = b
+c first solve ctrans(u)*y = b
+c
+ do 60 k = 1, n
+ t = zdotc(k-1,a(1,k),1,b(1),1)
+ b(k) = (b(k) - t)/dconjg(a(k,k))
+ 60 continue
+c
+c now solve ctrans(l)*x = y
+c
+ if (nm1 .lt. 1) go to 90
+ do 80 kb = 1, nm1
+ k = n - kb
+ b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1)
+ l = ipvt(k)
+ if (l .eq. k) go to 70
+ t = b(l)
+ b(l) = b(k)
+ b(k) = t
+ 70 continue
+ 80 continue
+ 90 continue
+ 100 continue
+ return
+ end
Property changes on: trunk/scipy/integrate/linpack_lite/zgesl.f
___________________________________________________________________
Name: svn:eol-style
+ native
Added: trunk/scipy/integrate/odepack/zvode.f
===================================================================
--- trunk/scipy/integrate/odepack/zvode.f 2008-02-20 05:27:16 UTC (rev 3950)
+++ trunk/scipy/integrate/odepack/zvode.f 2008-02-20 05:39:03 UTC (rev 3951)
@@ -0,0 +1,3650 @@
+*DECK ZVODE
+ SUBROUTINE ZVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
+ 1 ISTATE, IOPT, ZWORK, LZW, RWORK, LRW, IWORK, LIW,
+ 2 JAC, MF, RPAR, IPAR)
+ EXTERNAL F, JAC
+ DOUBLE COMPLEX Y, ZWORK
+ DOUBLE PRECISION T, TOUT, RTOL, ATOL, RWORK
+ INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LZW, LRW, IWORK, LIW,
+ 1 MF, IPAR
+ DIMENSION Y(*), RTOL(*), ATOL(*), ZWORK(LZW), RWORK(LRW),
+ 1 IWORK(LIW), RPAR(*), IPAR(*)
+C-----------------------------------------------------------------------
+C ZVODE: Variable-coefficient Ordinary Differential Equation solver,
+C with fixed-leading-coefficient implementation.
+C This version is in complex double precision.
+C
+C ZVODE solves the initial value problem for stiff or nonstiff
+C systems of first order ODEs,
+C dy/dt = f(t,y) , or, in component form,
+C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
+C Here the y vector is treated as complex.
+C ZVODE is a package based on the EPISODE and EPISODEB packages, and
+C on the ODEPACK user interface standard, with minor modifications.
+C
+C NOTE: When using ZVODE for a stiff system, it should only be used for
+C the case in which the function f is analytic, that is, when each f(i)
+C is an analytic function of each y(j). Analyticity means that the
+C partial derivative df(i)/dy(j) is a unique complex number, and this
+C fact is critical in the way ZVODE solves the dense or banded linear
+C systems that arise in the stiff case. For a complex stiff ODE system
+C in which f is not analytic, ZVODE is likely to have convergence
+C failures, and for this problem one should instead use DVODE on the
+C equivalent real system (in the real and imaginary parts of y).
+C-----------------------------------------------------------------------
+C Authors:
+C Peter N. Brown and Alan C. Hindmarsh
+C Center for Applied Scientific Computing
+C Lawrence Livermore National Laboratory
+C Livermore, CA 94551
+C and
+C George D. Byrne (Prof. Emeritus)
+C Illinois Institute of Technology
+C Chicago, IL 60616
+C-----------------------------------------------------------------------
+C For references, see DVODE.
+C-----------------------------------------------------------------------
+C Summary of usage.
+C
+C Communication between the user and the ZVODE package, for normal
+C situations, is summarized here. This summary describes only a subset
+C of the full set of options available. See the full description for
+C details, including optional communication, nonstandard options,
+C and instructions for special situations. See also the example
+C problem (with program and output) following this summary.
+C
+C A. First provide a subroutine of the form:
+C SUBROUTINE F (NEQ, T, Y, YDOT, RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), YDOT(NEQ)
+C DOUBLE PRECISION T
+C which supplies the vector function f by loading YDOT(i) with f(i).
+C
+C B. Next determine (or guess) whether or not the problem is stiff.
+C Stiffness occurs when the Jacobian matrix df/dy has an eigenvalue
+C whose real part is negative and large in magnitude, compared to the
+C reciprocal of the t span of interest. If the problem is nonstiff,
+C use a method flag MF = 10. If it is stiff, there are four standard
+C choices for MF (21, 22, 24, 25), and ZVODE requires the Jacobian
+C matrix in some form. In these cases (MF .gt. 0), ZVODE will use a
+C saved copy of the Jacobian matrix. If this is undesirable because of
+C storage limitations, set MF to the corresponding negative value
+C (-21, -22, -24, -25). (See full description of MF below.)
+C The Jacobian matrix is regarded either as full (MF = 21 or 22),
+C or banded (MF = 24 or 25). In the banded case, ZVODE requires two
+C half-bandwidth parameters ML and MU. These are, respectively, the
+C widths of the lower and upper parts of the band, excluding the main
+C diagonal. Thus the band consists of the locations (i,j) with
+C i-ML .le. j .le. i+MU, and the full bandwidth is ML+MU+1.
+C
+C C. If the problem is stiff, you are encouraged to supply the Jacobian
+C directly (MF = 21 or 24), but if this is not feasible, ZVODE will
+C compute it internally by difference quotients (MF = 22 or 25).
+C If you are supplying the Jacobian, provide a subroutine of the form:
+C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD, RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), PD(NROWPD,NEQ)
+C DOUBLE PRECISION T
+C which supplies df/dy by loading PD as follows:
+C For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j),
+C the partial derivative of f(i) with respect to y(j). (Ignore the
+C ML and MU arguments in this case.)
+C For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with
+C df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of
+C PD from the top down.
+C In either case, only nonzero elements need be loaded.
+C
+C D. Write a main program which calls subroutine ZVODE once for
+C each point at which answers are desired. This should also provide
+C for possible use of logical unit 6 for output of error messages
+C by ZVODE. On the first call to ZVODE, supply arguments as follows:
+C F = Name of subroutine for right-hand side vector f.
+C This name must be declared external in calling program.
+C NEQ = Number of first order ODEs.
+C Y = Double complex array of initial values, of length NEQ.
+C T = The initial value of the independent variable.
+C TOUT = First point where output is desired (.ne. T).
+C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
+C RTOL = Relative tolerance parameter (scalar).
+C ATOL = Absolute tolerance parameter (scalar or array).
+C The estimated local error in Y(i) will be controlled so as
+C to be roughly less (in magnitude) than
+C EWT(i) = RTOL*abs(Y(i)) + ATOL if ITOL = 1, or
+C EWT(i) = RTOL*abs(Y(i)) + ATOL(i) if ITOL = 2.
+C Thus the local error test passes if, in each component,
+C either the absolute error is less than ATOL (or ATOL(i)),
+C or the relative error is less than RTOL.
+C Use RTOL = 0.0 for pure absolute error control, and
+C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
+C control. Caution: Actual (global) errors may exceed these
+C local tolerances, so choose them conservatively.
+C ITASK = 1 for normal computation of output values of Y at t = TOUT.
+C ISTATE = Integer flag (input and output). Set ISTATE = 1.
+C IOPT = 0 to indicate no optional input used.
+C ZWORK = Double precision complex work array of length at least:
+C 15*NEQ for MF = 10,
+C 8*NEQ + 2*NEQ**2 for MF = 21 or 22,
+C 10*NEQ + (3*ML + 2*MU)*NEQ for MF = 24 or 25.
+C LZW = Declared length of ZWORK (in user's DIMENSION statement).
+C RWORK = Real work array of length at least 20 + NEQ.
+C LRW = Declared length of RWORK (in user's DIMENSION statement).
+C IWORK = Integer work array of length at least:
+C 30 for MF = 10,
+C 30 + NEQ for MF = 21, 22, 24, or 25.
+C If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower
+C and upper half-bandwidths ML,MU.
+C LIW = Declared length of IWORK (in user's DIMENSION statement).
+C JAC = Name of subroutine for Jacobian matrix (MF = 21 or 24).
+C If used, this name must be declared external in calling
+C program. If not used, pass a dummy name.
+C MF = Method flag. Standard values are:
+C 10 for nonstiff (Adams) method, no Jacobian used.
+C 21 for stiff (BDF) method, user-supplied full Jacobian.
+C 22 for stiff method, internally generated full Jacobian.
+C 24 for stiff method, user-supplied banded Jacobian.
+C 25 for stiff method, internally generated banded Jacobian.
+C RPAR = user-defined real or complex array passed to F and JAC.
+C IPAR = user-defined integer array passed to F and JAC.
+C Note that the main program must declare arrays Y, ZWORK, RWORK, IWORK,
+C and possibly ATOL, RPAR, and IPAR. RPAR may be declared REAL, DOUBLE,
+C COMPLEX, or DOUBLE COMPLEX, depending on the user's needs.
+C
+C E. The output from the first call (or any call) is:
+C Y = Array of computed values of y(t) vector.
+C T = Corresponding value of independent variable (normally TOUT).
+C ISTATE = 2 if ZVODE was successful, negative otherwise.
+C -1 means excess work done on this call. (Perhaps wrong MF.)
+C -2 means excess accuracy requested. (Tolerances too small.)
+C -3 means illegal input detected. (See printed message.)
+C -4 means repeated error test failures. (Check all input.)
+C -5 means repeated convergence failures. (Perhaps bad
+C Jacobian supplied or wrong choice of MF or tolerances.)
+C -6 means error weight became zero during problem. (Solution
+C component i vanished, and ATOL or ATOL(i) = 0.)
+C
+C F. To continue the integration after a successful return, simply
+C reset TOUT and call ZVODE again. No other parameters need be reset.
+C
+C-----------------------------------------------------------------------
+C EXAMPLE PROBLEM
+C
+C The program below uses ZVODE to solve the following system of 2 ODEs:
+C dw/dt = -i*w*w*z, dz/dt = i*z; w(0) = 1/2.1, z(0) = 1; t = 0 to 2*pi.
+C Solution: w = 1/(z + 1.1), z = exp(it). As z traces the unit circle,
+C w traces a circle of radius 10/2.1 with center at 11/2.1.
+C For convenience, Main passes RPAR = (imaginary unit i) to FEX and JEX.
+C
+C EXTERNAL FEX, JEX
+C DOUBLE COMPLEX Y(2), ZWORK(24), RPAR, WTRU, ERR
+C DOUBLE PRECISION ABERR, AEMAX, ATOL, RTOL, RWORK(22), T, TOUT
+C DIMENSION IWORK(32)
+C NEQ = 2
+C Y(1) = 1.0D0/2.1D0
+C Y(2) = 1.0D0
+C T = 0.0D0
+C DTOUT = 0.1570796326794896D0
+C TOUT = DTOUT
+C ITOL = 1
+C RTOL = 1.D-9
+C ATOL = 1.D-8
+C ITASK = 1
+C ISTATE = 1
+C IOPT = 0
+C LZW = 24
+C LRW = 22
+C LIW = 32
+C MF = 21
+C RPAR = DCMPLX(0.0D0,1.0D0)
+C AEMAX = 0.0D0
+C WRITE(6,10)
+C 10 FORMAT(' t',11X,'w',26X,'z')
+C DO 40 IOUT = 1,40
+C CALL ZVODE(FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,
+C 1 ZWORK,LZW,RWORK,LRW,IWORK,LIW,JEX,MF,RPAR,IPAR)
+C WTRU = 1.0D0/DCMPLX(COS(T) + 1.1D0, SIN(T))
+C ERR = Y(1) - WTRU
+C ABERR = ABS(DREAL(ERR)) + ABS(DIMAG(ERR))
+C AEMAX = MAX(AEMAX,ABERR)
+C WRITE(6,20) T, DREAL(Y(1)),DIMAG(Y(1)), DREAL(Y(2)),DIMAG(Y(2))
+C 20 FORMAT(F9.5,2X,2F12.7,3X,2F12.7)
+C IF (ISTATE .LT. 0) THEN
+C WRITE(6,30) ISTATE
+C 30 FORMAT(//'***** Error halt. ISTATE =',I3)
+C STOP
+C ENDIF
+C 40 TOUT = TOUT + DTOUT
+C WRITE(6,50) IWORK(11), IWORK(12), IWORK(13), IWORK(20),
+C 1 IWORK(21), IWORK(22), IWORK(23), AEMAX
+C 50 FORMAT(/' No. steps =',I4,' No. f-s =',I5,
+C 1 ' No. J-s =',I4,' No. LU-s =',I4/
+C 2 ' No. nonlinear iterations =',I4/
+C 3 ' No. nonlinear convergence failures =',I4/
+C 4 ' No. error test failures =',I4/
+C 5 ' Max. abs. error in w =',D10.2)
+C STOP
+C END
+C
+C SUBROUTINE FEX (NEQ, T, Y, YDOT, RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), YDOT(NEQ), RPAR
+C DOUBLE PRECISION T
+C YDOT(1) = -RPAR*Y(1)*Y(1)*Y(2)
+C YDOT(2) = RPAR*Y(2)
+C RETURN
+C END
+C
+C SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD, RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), PD(NRPD,NEQ), RPAR
+C DOUBLE PRECISION T
+C PD(1,1) = -2.0D0*RPAR*Y(1)*Y(2)
+C PD(1,2) = -RPAR*Y(1)*Y(1)
+C PD(2,2) = RPAR
+C RETURN
+C END
+C
+C The output of this example program is as follows:
+C
+C t w z
+C 0.15708 0.4763242 -0.0356919 0.9876884 0.1564345
+C 0.31416 0.4767322 -0.0718256 0.9510565 0.3090170
+C 0.47124 0.4774351 -0.1088651 0.8910065 0.4539906
+C 0.62832 0.4784699 -0.1473206 0.8090170 0.5877853
+C 0.78540 0.4798943 -0.1877789 0.7071067 0.7071069
+C 0.94248 0.4817938 -0.2309414 0.5877852 0.8090171
+C 1.09956 0.4842934 -0.2776778 0.4539904 0.8910066
+C 1.25664 0.4875766 -0.3291039 0.3090169 0.9510566
+C 1.41372 0.4919177 -0.3866987 0.1564343 0.9876884
+C 1.57080 0.4977376 -0.4524889 -0.0000001 1.0000000
+C 1.72788 0.5057044 -0.5293524 -0.1564346 0.9876883
+C 1.88496 0.5169274 -0.6215400 -0.3090171 0.9510565
+C 2.04204 0.5333540 -0.7356275 -0.4539906 0.8910065
+C 2.19911 0.5586542 -0.8823669 -0.5877854 0.8090169
+C 2.35619 0.6004188 -1.0806013 -0.7071069 0.7071067
+C 2.51327 0.6764486 -1.3664281 -0.8090171 0.5877851
+C 2.67035 0.8366909 -1.8175245 -0.8910066 0.4539904
+C 2.82743 1.2657121 -2.6260146 -0.9510566 0.3090168
+C 2.98451 3.0284506 -4.2182180 -0.9876884 0.1564343
+C 3.14159 10.0000699 0.0000663 -1.0000000 -0.0000002
+C 3.29867 3.0284170 4.2182053 -0.9876883 -0.1564346
+C 3.45575 1.2657041 2.6260067 -0.9510565 -0.3090172
+C 3.61283 0.8366878 1.8175205 -0.8910064 -0.4539907
+C 3.76991 0.6764469 1.3664259 -0.8090169 -0.5877854
+C 3.92699 0.6004178 1.0806000 -0.7071066 -0.7071069
+C 4.08407 0.5586535 0.8823662 -0.5877851 -0.8090171
+C 4.24115 0.5333535 0.7356271 -0.4539903 -0.8910066
+C 4.39823 0.5169271 0.6215398 -0.3090168 -0.9510566
+C 4.55531 0.5057041 0.5293523 -0.1564343 -0.9876884
+C 4.71239 0.4977374 0.4524890 0.0000002 -1.0000000
+C 4.86947 0.4919176 0.3866988 0.1564347 -0.9876883
+C 5.02655 0.4875765 0.3291040 0.3090172 -0.9510564
+C 5.18363 0.4842934 0.2776780 0.4539907 -0.8910064
+C 5.34071 0.4817939 0.2309415 0.5877854 -0.8090169
+C 5.49779 0.4798944 0.1877791 0.7071069 -0.7071066
+C 5.65487 0.4784700 0.1473208 0.8090171 -0.5877850
+C 5.81195 0.4774352 0.1088652 0.8910066 -0.4539903
+C 5.96903 0.4767324 0.0718257 0.9510566 -0.3090168
+C 6.12611 0.4763244 0.0356920 0.9876884 -0.1564342
+C 6.28319 0.4761907 0.0000000 1.0000000 0.0000003
+C
+C No. steps = 542 No. f-s = 610 No. J-s = 10 No. LU-s = 47
+C No. nonlinear iterations = 607
+C No. nonlinear convergence failures = 0
+C No. error test failures = 13
+C Max. abs. error in w = 0.13E-03
+C
+C-----------------------------------------------------------------------
+C Full description of user interface to ZVODE.
+C
+C The user interface to ZVODE consists of the following parts.
+C
+C i. The call sequence to subroutine ZVODE, which is a driver
+C routine for the solver. This includes descriptions of both
+C the call sequence arguments and of user-supplied routines.
+C Following these descriptions is
+C * a description of optional input available through the
+C call sequence,
+C * a description of optional output (in the work arrays), and
+C * instructions for interrupting and restarting a solution.
+C
+C ii. Descriptions of other routines in the ZVODE package that may be
+C (optionally) called by the user. These provide the ability to
+C alter error message handling, save and restore the internal
+C COMMON, and obtain specified derivatives of the solution y(t).
+C
+C iii. Descriptions of COMMON blocks to be declared in overlay
+C or similar environments.
+C
+C iv. Description of two routines in the ZVODE package, either of
+C which the user may replace with his own version, if desired.
+C these relate to the measurement of errors.
+C
+C-----------------------------------------------------------------------
+C Part i. Call Sequence.
+C
+C The call sequence parameters used for input only are
+C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
+C and those used for both input and output are
+C Y, T, ISTATE.
+C The work arrays ZWORK, RWORK, and IWORK are also used for conditional
+C and optional input and optional output. (The term output here refers
+C to the return from subroutine ZVODE to the user's calling program.)
+C
+C The legality of input parameters will be thoroughly checked on the
+C initial call for the problem, but not checked thereafter unless a
+C change in input parameters is flagged by ISTATE = 3 in the input.
+C
+C The descriptions of the call arguments are as follows.
+C
+C F = The name of the user-supplied subroutine defining the
+C ODE system. The system must be put in the first-order
+C form dy/dt = f(t,y), where f is a vector-valued function
+C of the scalar t and the vector y. Subroutine F is to
+C compute the function f. It is to have the form
+C SUBROUTINE F (NEQ, T, Y, YDOT, RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), YDOT(NEQ)
+C DOUBLE PRECISION T
+C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
+C is output. Y and YDOT are double complex arrays of length
+C NEQ. Subroutine F should not alter Y(1),...,Y(NEQ).
+C F must be declared EXTERNAL in the calling program.
+C
+C Subroutine F may access user-defined real/complex and
+C integer work arrays RPAR and IPAR, which are to be
+C dimensioned in the calling program.
+C
+C If quantities computed in the F routine are needed
+C externally to ZVODE, an extra call to F should be made
+C for this purpose, for consistent and accurate results.
+C If only the derivative dy/dt is needed, use ZVINDY instead.
+C
+C NEQ = The size of the ODE system (number of first order
+C ordinary differential equations). Used only for input.
+C NEQ may not be increased during the problem, but
+C can be decreased (with ISTATE = 3 in the input).
+C
+C Y = A double precision complex array for the vector of dependent
+C variables, of length NEQ or more. Used for both input and
+C output on the first call (ISTATE = 1), and only for output
+C on other calls. On the first call, Y must contain the
+C vector of initial values. In the output, Y contains the
+C computed solution evaluated at T. If desired, the Y array
+C may be used for other purposes between calls to the solver.
+C
+C This array is passed as the Y argument in all calls to
+C F and JAC.
+C
+C T = The independent variable. In the input, T is used only on
+C the first call, as the initial point of the integration.
+C In the output, after each call, T is the value at which a
+C computed solution Y is evaluated (usually the same as TOUT).
+C On an error return, T is the farthest point reached.
+C
+C TOUT = The next value of t at which a computed solution is desired.
+C Used only for input.
+C
+C When starting the problem (ISTATE = 1), TOUT may be equal
+C to T for one call, then should .ne. T for the next call.
+C For the initial T, an input value of TOUT .ne. T is used
+C in order to determine the direction of the integration
+C (i.e. the algebraic sign of the step sizes) and the rough
+C scale of the problem. Integration in either direction
+C (forward or backward in t) is permitted.
+C
+C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
+C the first call (i.e. the first call with TOUT .ne. T).
+C Otherwise, TOUT is required on every call.
+C
+C If ITASK = 1, 3, or 4, the values of TOUT need not be
+C monotone, but a value of TOUT which backs up is limited
+C to the current internal t interval, whose endpoints are
+C TCUR - HU and TCUR. (See optional output, below, for
+C TCUR and HU.)
+C
+C ITOL = An indicator for the type of error control. See
+C description below under ATOL. Used only for input.
+C
+C RTOL = A relative error tolerance parameter, either a scalar or
+C an array of length NEQ. See description below under ATOL.
+C Input only.
+C
+C ATOL = An absolute error tolerance parameter, either a scalar or
+C an array of length NEQ. Input only.
+C
+C The input parameters ITOL, RTOL, and ATOL determine
+C the error control performed by the solver. The solver will
+C control the vector e = (e(i)) of estimated local errors
+C in Y, according to an inequality of the form
+C rms-norm of ( e(i)/EWT(i) ) .le. 1,
+C where EWT(i) = RTOL(i)*abs(Y(i)) + ATOL(i),
+C and the rms-norm (root-mean-square norm) here is
+C rms-norm(v) = sqrt(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
+C is a vector of weights which must always be positive, and
+C the values of RTOL and ATOL should all be non-negative.
+C The following table gives the types (scalar/array) of
+C RTOL and ATOL, and the corresponding form of EWT(i).
+C
+C ITOL RTOL ATOL EWT(i)
+C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
+C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
+C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
+C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i)
+C
+C When either of these parameters is a scalar, it need not
+C be dimensioned in the user's calling program.
+C
+C If none of the above choices (with ITOL, RTOL, and ATOL
+C fixed throughout the problem) is suitable, more general
+C error controls can be obtained by substituting
+C user-supplied routines for the setting of EWT and/or for
+C the norm calculation. See Part iv below.
+C
+C If global errors are to be estimated by making a repeated
+C run on the same problem with smaller tolerances, then all
+C components of RTOL and ATOL (i.e. of EWT) should be scaled
+C down uniformly.
+C
+C ITASK = An index specifying the task to be performed.
+C Input only. ITASK has the following values and meanings.
+C 1 means normal computation of output values of y(t) at
+C t = TOUT (by overshooting and interpolating).
+C 2 means take one step only and return.
+C 3 means stop at the first internal mesh point at or
+C beyond t = TOUT and return.
+C 4 means normal computation of output values of y(t) at
+C t = TOUT but without overshooting t = TCRIT.
+C TCRIT must be input as RWORK(1). TCRIT may be equal to
+C or beyond TOUT, but not behind it in the direction of
+C integration. This option is useful if the problem
+C has a singularity at or beyond t = TCRIT.
+C 5 means take one step, without passing TCRIT, and return.
+C TCRIT must be input as RWORK(1).
+C
+C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
+C (within roundoff), it will return T = TCRIT (exactly) to
+C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
+C in which case answers at T = TOUT are returned first).
+C
+C ISTATE = an index used for input and output to specify the
+C the state of the calculation.
+C
+C In the input, the values of ISTATE are as follows.
+C 1 means this is the first call for the problem
+C (initializations will be done). See note below.
+C 2 means this is not the first call, and the calculation
+C is to continue normally, with no change in any input
+C parameters except possibly TOUT and ITASK.
+C (If ITOL, RTOL, and/or ATOL are changed between calls
+C with ISTATE = 2, the new values will be used but not
+C tested for legality.)
+C 3 means this is not the first call, and the
+C calculation is to continue normally, but with
+C a change in input parameters other than
+C TOUT and ITASK. Changes are allowed in
+C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU,
+C and any of the optional input except H0.
+C (See IWORK description for ML and MU.)
+C Note: A preliminary call with TOUT = T is not counted
+C as a first call here, as no initialization or checking of
+C input is done. (Such a call is sometimes useful to include
+C the initial conditions in the output.)
+C Thus the first call for which TOUT .ne. T requires
+C ISTATE = 1 in the input.
+C
+C In the output, ISTATE has the following values and meanings.
+C 1 means nothing was done, as TOUT was equal to T with
+C ISTATE = 1 in the input.
+C 2 means the integration was performed successfully.
+C -1 means an excessive amount of work (more than MXSTEP
+C steps) was done on this call, before completing the
+C requested task, but the integration was otherwise
+C successful as far as T. (MXSTEP is an optional input
+C and is normally 500.) To continue, the user may
+C simply reset ISTATE to a value .gt. 1 and call again.
+C (The excess work step counter will be reset to 0.)
+C In addition, the user may increase MXSTEP to avoid
+C this error return. (See optional input below.)
+C -2 means too much accuracy was requested for the precision
+C of the machine being used. This was detected before
+C completing the requested task, but the integration
+C was successful as far as T. To continue, the tolerance
+C parameters must be reset, and ISTATE must be set
+C to 3. The optional output TOLSF may be used for this
+C purpose. (Note: If this condition is detected before
+C taking any steps, then an illegal input return
+C (ISTATE = -3) occurs instead.)
+C -3 means illegal input was detected, before taking any
+C integration steps. See written message for details.
+C Note: If the solver detects an infinite loop of calls
+C to the solver with illegal input, it will cause
+C the run to stop.
+C -4 means there were repeated error test failures on
+C one attempted step, before completing the requested
+C task, but the integration was successful as far as T.
+C The problem may have a singularity, or the input
+C may be inappropriate.
+C -5 means there were repeated convergence test failures on
+C one attempted step, before completing the requested
+C task, but the integration was successful as far as T.
+C This may be caused by an inaccurate Jacobian matrix,
+C if one is being used.
+C -6 means EWT(i) became zero for some i during the
+C integration. Pure relative error control (ATOL(i)=0.0)
+C was requested on a variable which has now vanished.
+C The integration was successful as far as T.
+C
+C Note: Since the normal output value of ISTATE is 2,
+C it does not need to be reset for normal continuation.
+C Also, since a negative input value of ISTATE will be
+C regarded as illegal, a negative output value requires the
+C user to change it, and possibly other input, before
+C calling the solver again.
+C
+C IOPT = An integer flag to specify whether or not any optional
+C input is being used on this call. Input only.
+C The optional input is listed separately below.
+C IOPT = 0 means no optional input is being used.
+C Default values will be used in all cases.
+C IOPT = 1 means optional input is being used.
+C
+C ZWORK = A double precision complex working array.
+C The length of ZWORK must be at least
+C NYH*(MAXORD + 1) + 2*NEQ + LWM where
+C NYH = the initial value of NEQ,
+C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
+C smaller value is given as an optional input),
+C LWM = length of work space for matrix-related data:
+C LWM = 0 if MITER = 0,
+C LWM = 2*NEQ**2 if MITER = 1 or 2, and MF.gt.0,
+C LWM = NEQ**2 if MITER = 1 or 2, and MF.lt.0,
+C LWM = NEQ if MITER = 3,
+C LWM = (3*ML+2*MU+2)*NEQ if MITER = 4 or 5, and MF.gt.0,
+C LWM = (2*ML+MU+1)*NEQ if MITER = 4 or 5, and MF.lt.0.
+C (See the MF description for METH and MITER.)
+C Thus if MAXORD has its default value and NEQ is constant,
+C this length is:
+C 15*NEQ for MF = 10,
+C 15*NEQ + 2*NEQ**2 for MF = 11 or 12,
+C 15*NEQ + NEQ**2 for MF = -11 or -12,
+C 16*NEQ for MF = 13,
+C 17*NEQ + (3*ML+2*MU)*NEQ for MF = 14 or 15,
+C 16*NEQ + (2*ML+MU)*NEQ for MF = -14 or -15,
+C 8*NEQ for MF = 20,
+C 8*NEQ + 2*NEQ**2 for MF = 21 or 22,
+C 8*NEQ + NEQ**2 for MF = -21 or -22,
+C 9*NEQ for MF = 23,
+C 10*NEQ + (3*ML+2*MU)*NEQ for MF = 24 or 25.
+C 9*NEQ + (2*ML+MU)*NEQ for MF = -24 or -25.
+C
+C LZW = The length of the array ZWORK, as declared by the user.
+C (This will be checked by the solver.)
+C
+C RWORK = A real working array (double precision).
+C The length of RWORK must be at least 20 + NEQ.
+C The first 20 words of RWORK are reserved for conditional
+C and optional input and optional output.
+C
+C The following word in RWORK is a conditional input:
+C RWORK(1) = TCRIT = critical value of t which the solver
+C is not to overshoot. Required if ITASK is
+C 4 or 5, and ignored otherwise. (See ITASK.)
+C
+C LRW = The length of the array RWORK, as declared by the user.
+C (This will be checked by the solver.)
+C
+C IWORK = An integer work array. The length of IWORK must be at least
+C 30 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or
+C 30 + NEQ otherwise (abs(MF) = 11,12,14,15,21,22,24,25).
+C The first 30 words of IWORK are reserved for conditional and
+C optional input and optional output.
+C
+C The following 2 words in IWORK are conditional input:
+C IWORK(1) = ML These are the lower and upper
+C IWORK(2) = MU half-bandwidths, respectively, of the
+C banded Jacobian, excluding the main diagonal.
+C The band is defined by the matrix locations
+C (i,j) with i-ML .le. j .le. i+MU. ML and MU
+C must satisfy 0 .le. ML,MU .le. NEQ-1.
+C These are required if MITER is 4 or 5, and
+C ignored otherwise. ML and MU may in fact be
+C the band parameters for a matrix to which
+C df/dy is only approximately equal.
+C
+C LIW = the length of the array IWORK, as declared by the user.
+C (This will be checked by the solver.)
+C
+C Note: The work arrays must not be altered between calls to ZVODE
+C for the same problem, except possibly for the conditional and
+C optional input, and except for the last 2*NEQ words of ZWORK and
+C the last NEQ words of RWORK. The latter space is used for internal
+C scratch space, and so is available for use by the user outside ZVODE
+C between calls, if desired (but not for use by F or JAC).
+C
+C JAC = The name of the user-supplied routine (MITER = 1 or 4) to
+C compute the Jacobian matrix, df/dy, as a function of
+C the scalar t and the vector y. It is to have the form
+C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD,
+C RPAR, IPAR)
+C DOUBLE COMPLEX Y(NEQ), PD(NROWPD,NEQ)
+C DOUBLE PRECISION T
+C where NEQ, T, Y, ML, MU, and NROWPD are input and the array
+C PD is to be loaded with partial derivatives (elements of the
+C Jacobian matrix) in the output. PD must be given a first
+C dimension of NROWPD. T and Y have the same meaning as in
+C Subroutine F.
+C In the full matrix case (MITER = 1), ML and MU are
+C ignored, and the Jacobian is to be loaded into PD in
+C columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
+C In the band matrix case (MITER = 4), the elements
+C within the band are to be loaded into PD in columnwise
+C manner, with diagonal lines of df/dy loaded into the rows
+C of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j).
+C ML and MU are the half-bandwidth parameters. (See IWORK).
+C The locations in PD in the two triangular areas which
+C correspond to nonexistent matrix elements can be ignored
+C or loaded arbitrarily, as they are overwritten by ZVODE.
+C JAC need not provide df/dy exactly. A crude
+C approximation (possibly with a smaller bandwidth) will do.
+C In either case, PD is preset to zero by the solver,
+C so that only the nonzero elements need be loaded by JAC.
+C Each call to JAC is preceded by a call to F with the same
+C arguments NEQ, T, and Y. Thus to gain some efficiency,
+C intermediate quantities shared by both calculations may be
+C saved in a user COMMON block by F and not recomputed by JAC,
+C if desired. Also, JAC may alter the Y array, if desired.
+C JAC must be declared external in the calling program.
+C Subroutine JAC may access user-defined real/complex and
+C integer work arrays, RPAR and IPAR, whose dimensions are set
+C by the user in the calling program.
+C
+C MF = The method flag. Used only for input. The legal values of
+C MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25,
+C -11, -12, -14, -15, -21, -22, -24, -25.
+C MF is a signed two-digit integer, MF = JSV*(10*METH + MITER).
+C JSV = SIGN(MF) indicates the Jacobian-saving strategy:
+C JSV = 1 means a copy of the Jacobian is saved for reuse
+C in the corrector iteration algorithm.
+C JSV = -1 means a copy of the Jacobian is not saved
+C (valid only for MITER = 1, 2, 4, or 5).
+C METH indicates the basic linear multistep method:
+C METH = 1 means the implicit Adams method.
+C METH = 2 means the method based on backward
+C differentiation formulas (BDF-s).
+C MITER indicates the corrector iteration method:
+C MITER = 0 means functional iteration (no Jacobian matrix
+C is involved).
+C MITER = 1 means chord iteration with a user-supplied
+C full (NEQ by NEQ) Jacobian.
+C MITER = 2 means chord iteration with an internally
+C generated (difference quotient) full Jacobian
+C (using NEQ extra calls to F per df/dy value).
+C MITER = 3 means chord iteration with an internally
+C generated diagonal Jacobian approximation
+C (using 1 extra call to F per df/dy evaluation).
+C MITER = 4 means chord iteration with a user-supplied
+C banded Jacobian.
+C MITER = 5 means chord iteration with an internally
+C generated banded Jacobian (using ML+MU+1 extra
+C calls to F per df/dy evaluation).
+C If MITER = 1 or 4, the user must supply a subroutine JAC
+C (the name is arbitrary) as described above under JAC.
+C For other values of MITER, a dummy argument can be used.
+C
+C RPAR User-specified array used to communicate real or complex
+C parameters to user-supplied subroutines. If RPAR is an
+C array, it must be dimensioned in the user's calling program;
+C if it is unused or it is a scalar, then it need not be
+C dimensioned. The type of RPAR may be REAL, DOUBLE, COMPLEX,
+C or DOUBLE COMPLEX, depending on the user program's needs.
+C RPAR is not type-declared within ZVODE, but simply passed
+C (by address) to the user's F and JAC routines.
+C
+C IPAR User-specified array used to communicate integer parameter
+C to user-supplied subroutines. If IPAR is an array, it must
+C be dimensioned in the user's calling program.
+C-----------------------------------------------------------------------
+C Optional Input.
+C
+C The following is a list of the optional input provided for in the
+C call sequence. (See also Part ii.) For each such input variable,
+C this table lists its name as used in this documentation, its
+C location in the call sequence, its meaning, and the default value.
+C The use of any of this input requires IOPT = 1, and in that
+C case all of this input is examined. A value of zero for any
+C of these optional input variables will cause the default value to be
+C used. Thus to use a subset of the optional input, simply preload
+C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
+C then set those of interest to nonzero values.
+C
+C NAME LOCATION MEANING AND DEFAULT VALUE
+C
+C H0 RWORK(5) The step size to be attempted on the first step.
+C The default value is determined by the solver.
+C
+C HMAX RWORK(6) The maximum absolute step size allowed.
+C The default value is infinite.
+C
+C HMIN RWORK(7) The minimum absolute step size allowed.
+C The default value is 0. (This lower bound is not
+C enforced on the final step before reaching TCRIT
+C when ITASK = 4 or 5.)
+C
+C MAXORD IWORK(5) The maximum order to be allowed. The default
+C value is 12 if METH = 1, and 5 if METH = 2.
+C If MAXORD exceeds the default value, it will
+C be reduced to the default value.
+C If MAXORD is changed during the problem, it may
+C cause the current order to be reduced.
+C
+C MXSTEP IWORK(6) Maximum number of (internally defined) steps
+C allowed during one call to the solver.
+C The default value is 500.
+C
+C MXHNIL IWORK(7) Maximum number of messages printed (per problem)
+C warning that T + H = T on a step (H = step size).
+C This must be positive to result in a non-default
+C value. The default value is 10.
+C
+C-----------------------------------------------------------------------
+C Optional Output.
+C
+C As optional additional output from ZVODE, the variables listed
+C below are quantities related to the performance of ZVODE
+C which are available to the user. These are communicated by way of
+C the work arrays, but also have internal mnemonic names as shown.
+C Except where stated otherwise, all of this output is defined
+C on any successful return from ZVODE, and on any return with
+C ISTATE = -1, -2, -4, -5, or -6. On an illegal input return
+C (ISTATE = -3), they will be unchanged from their existing values
+C (if any), except possibly for TOLSF, LENZW, LENRW, and LENIW.
+C On any error return, output relevant to the error will be defined,
+C as noted below.
+C
+C NAME LOCATION MEANING
+C
+C HU RWORK(11) The step size in t last used (successfully).
+C
+C HCUR RWORK(12) The step size to be attempted on the next step.
+C
+C TCUR RWORK(13) The current value of the independent variable
+C which the solver has actually reached, i.e. the
+C current internal mesh point in t. In the output,
+C TCUR will always be at least as far from the
+C initial value of t as the current argument T,
+C but may be farther (if interpolation was done).
+C
+C TOLSF RWORK(14) A tolerance scale factor, greater than 1.0,
+C computed when a request for too much accuracy was
+C detected (ISTATE = -3 if detected at the start of
+C the problem, ISTATE = -2 otherwise). If ITOL is
+C left unaltered but RTOL and ATOL are uniformly
+C scaled up by a factor of TOLSF for the next call,
+C then the solver is deemed likely to succeed.
+C (The user may also ignore TOLSF and alter the
+C tolerance parameters in any other way appropriate.)
+C
+C NST IWORK(11) The number of steps taken for the problem so far.
+C
+C NFE IWORK(12) The number of f evaluations for the problem so far.
+C
+C NJE IWORK(13) The number of Jacobian evaluations so far.
+C
+C NQU IWORK(14) The method order last used (successfully).
+C
+C NQCUR IWORK(15) The order to be attempted on the next step.
+C
+C IMXER IWORK(16) The index of the component of largest magnitude in
+C the weighted local error vector ( e(i)/EWT(i) ),
+C on an error return with ISTATE = -4 or -5.
+C
+C LENZW IWORK(17) The length of ZWORK actually required.
+C This is defined on normal returns and on an illegal
+C input return for insufficient storage.
+C
+C LENRW IWORK(18) The length of RWORK actually required.
+C This is defined on normal returns and on an illegal
+C input return for insufficient storage.
+C
+C LENIW IWORK(19) The length of IWORK actually required.
+C This is defined on normal returns and on an illegal
+C input return for insufficient storage.
+C
+C NLU IWORK(20) The number of matrix LU decompositions so far.
+C
+C NNI IWORK(21) The number of nonlinear (Newton) iterations so far.
+C
+C NCFN IWORK(22) The number of convergence failures of the nonlinear
+C solver so far.
+C
+C NETF IWORK(23) The number of error test failures of the integrator
+C so far.
+C
+C The following two arrays are segments of the ZWORK array which
+C may also be of interest to the user as optional output.
+C For each array, the table below gives its internal name,
+C its base address in ZWORK, and its description.
+C
+C NAME BASE ADDRESS DESCRIPTION
+C
+C YH 1 The Nordsieck history array, of size NYH by
+C (NQCUR + 1), where NYH is the initial value
+C of NEQ. For j = 0,1,...,NQCUR, column j+1
+C of YH contains HCUR**j/factorial(j) times
+C the j-th derivative of the interpolating
+C polynomial currently representing the
+C solution, evaluated at t = TCUR.
+C
+C ACOR LENZW-NEQ+1 Array of size NEQ used for the accumulated
+C corrections on each step, scaled in the output
+C to represent the estimated local error in Y
+C on the last step. This is the vector e in
+C the description of the error control. It is
+C defined only on a successful return from ZVODE.
+C
+C-----------------------------------------------------------------------
+C Interrupting and Restarting
+C
+C If the integration of a given problem by ZVODE is to be
+C interrrupted and then later continued, such as when restarting
+C an interrupted run or alternating between two or more ODE problems,
+C the user should save, following the return from the last ZVODE call
+C prior to the interruption, the contents of the call sequence
+C variables and internal COMMON blocks, and later restore these
+C values before the next ZVODE call for that problem. To save
+C and restore the COMMON blocks, use subroutine ZVSRCO, as
+C described below in part ii.
+C
+C In addition, if non-default values for either LUN or MFLAG are
+C desired, an extra call to XSETUN and/or XSETF should be made just
+C before continuing the integration. See Part ii below for details.
+C
+C-----------------------------------------------------------------------
+C Part ii. Other Routines Callable.
+C
+C The following are optional calls which the user may make to
+C gain additional capabilities in conjunction with ZVODE.
+C (The routines XSETUN and XSETF are designed to conform to the
+C SLATEC error handling package.)
+C
+C FORM OF CALL FUNCTION
+C CALL XSETUN(LUN) Set the logical unit number, LUN, for
+C output of messages from ZVODE, if
+C the default is not desired.
+C The default value of LUN is 6.
+C
+C CALL XSETF(MFLAG) Set a flag to control the printing of
+C messages by ZVODE.
+C MFLAG = 0 means do not print. (Danger:
+C This risks losing valuable information.)
+C MFLAG = 1 means print (the default).
+C
+C Either of the above calls may be made at
+C any time and will take effect immediately.
+C
+C CALL ZVSRCO(RSAV,ISAV,JOB) Saves and restores the contents of
+C the internal COMMON blocks used by
+C ZVODE. (See Part iii below.)
+C RSAV must be a real array of length 51
+C or more, and ISAV must be an integer
+C array of length 40 or more.
+C JOB=1 means save COMMON into RSAV/ISAV.
+C JOB=2 means restore COMMON from RSAV/ISAV.
+C ZVSRCO is useful if one is
+C interrupting a run and restarting
+C later, or alternating between two or
+C more problems solved with ZVODE.
+C
+C CALL ZVINDY(,,,,,) Provide derivatives of y, of various
+C (See below.) orders, at a specified point T, if
+C desired. It may be called only after
+C a successful return from ZVODE.
+C
+C The detailed instructions for using ZVINDY are as follows.
+C The form of the call is:
+C
+C CALL ZVINDY (T, K, ZWORK, NYH, DKY, IFLAG)
+C
+C The input parameters are:
+C
+C T = Value of independent variable where answers are desired
+C (normally the same as the T last returned by ZVODE).
+C For valid results, T must lie between TCUR - HU and TCUR.
+C (See optional output for TCUR and HU.)
+C K = Integer order of the derivative desired. K must satisfy
+C 0 .le. K .le. NQCUR, where NQCUR is the current order
+C (see optional output). The capability corresponding
+C to K = 0, i.e. computing y(T), is already provided
+C by ZVODE directly. Since NQCUR .ge. 1, the first
+C derivative dy/dt is always available with ZVINDY.
+C ZWORK = The history array YH.
+C NYH = Column length of YH, equal to the initial value of NEQ.
+C
+C The output parameters are:
+C
+C DKY = A double complex array of length NEQ containing the
+C computed value of the K-th derivative of y(t).
+C IFLAG = Integer flag, returned as 0 if K and T were legal,
+C -1 if K was illegal, and -2 if T was illegal.
+C On an error return, a message is also written.
+C-----------------------------------------------------------------------
+C Part iii. COMMON Blocks.
+C If ZVODE is to be used in an overlay situation, the user
+C must declare, in the primary overlay, the variables in:
+C (1) the call sequence to ZVODE,
+C (2) the two internal COMMON blocks
+C /ZVOD01/ of length 83 (50 double precision words
+C followed by 33 integer words),
+C /ZVOD02/ of length 9 (1 double precision word
+C followed by 8 integer words),
+C
+C If ZVODE is used on a system in which the contents of internal
+C COMMON blocks are not preserved between calls, the user should
+C declare the above two COMMON blocks in his calling program to insure
+C that their contents are preserved.
+C
+C-----------------------------------------------------------------------
+C Part iv. Optionally Replaceable Solver Routines.
+C
+C Below are descriptions of two routines in the ZVODE package which
+C relate to the measurement of errors. Either routine can be
+C replaced by a user-supplied version, if desired. However, since such
+C a replacement may have a major impact on performance, it should be
+C done only when absolutely necessary, and only with great caution.
+C (Note: The means by which the package version of a routine is
+C superseded by the user's version may be system-dependent.)
+C
+C (a) ZEWSET.
+C The following subroutine is called just before each internal
+C integration step, and sets the array of error weights, EWT, as
+C described under ITOL/RTOL/ATOL above:
+C SUBROUTINE ZEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
+C where NEQ, ITOL, RTOL, and ATOL are as in the ZVODE call sequence,
+C YCUR contains the current (double complex) dependent variable vector,
+C and EWT is the array of weights set by ZEWSET.
+C
+C If the user supplies this subroutine, it must return in EWT(i)
+C (i = 1,...,NEQ) a positive quantity suitable for comparison with
+C errors in Y(i). The EWT array returned by ZEWSET is passed to the
+C ZVNORM routine (See below.), and also used by ZVODE in the computation
+C of the optional output IMXER, the diagonal Jacobian approximation,
+C and the increments for difference quotient Jacobians.
+C
+C In the user-supplied version of ZEWSET, it may be desirable to use
+C the current values of derivatives of y. Derivatives up to order NQ
+C are available from the history array YH, described above under
+C Optional Output. In ZEWSET, YH is identical to the YCUR array,
+C extended to NQ + 1 columns with a column length of NYH and scale
+C factors of h**j/factorial(j). On the first call for the problem,
+C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
+C NYH is the initial value of NEQ. The quantities NQ, H, and NST
+C can be obtained by including in ZEWSET the statements:
+C DOUBLE PRECISION RVOD, H, HU
+C COMMON /ZVOD01/ RVOD(50), IVOD(33)
+C COMMON /ZVOD02/ HU, NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C NQ = IVOD(28)
+C H = RVOD(21)
+C Thus, for example, the current value of dy/dt can be obtained as
+C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
+C unnecessary when NST = 0).
+C
+C (b) ZVNORM.
+C The following is a real function routine which computes the weighted
+C root-mean-square norm of a vector v:
+C D = ZVNORM (N, V, W)
+C where:
+C N = the length of the vector,
+C V = double complex array of length N containing the vector,
+C W = real array of length N containing weights,
+C D = sqrt( (1/N) * sum(abs(V(i))*W(i))**2 ).
+C ZVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
+C EWT is as set by subroutine ZEWSET.
+C
+C If the user supplies this function, it should return a non-negative
+C value of ZVNORM suitable for use in the error control in ZVODE.
+C None of the arguments should be altered by ZVNORM.
+C For example, a user-supplied ZVNORM routine might:
+C -substitute a max-norm of (V(i)*W(i)) for the rms-norm, or
+C -ignore some components of V in the norm, with the effect of
+C suppressing the error control on those components of Y.
+C-----------------------------------------------------------------------
+C REVISION HISTORY (YYYYMMDD)
+C 20060517 DATE WRITTEN, modified from DVODE of 20020430.
+C 20061227 Added note on use for analytic f.
+C-----------------------------------------------------------------------
+C Other Routines in the ZVODE Package.
+C
+C In addition to Subroutine ZVODE, the ZVODE package includes the
+C following subroutines and function routines:
+C ZVHIN computes an approximate step size for the initial step.
+C ZVINDY computes an interpolated value of the y vector at t = TOUT.
+C ZVSTEP is the core integrator, which does one step of the
+C integration and the associated error control.
+C ZVSET sets all method coefficients and test constants.
+C ZVNLSD solves the underlying nonlinear system -- the corrector.
+C ZVJAC computes and preprocesses the Jacobian matrix J = df/dy
+C and the Newton iteration matrix P = I - (h/l1)*J.
+C ZVSOL manages solution of linear system in chord iteration.
+C ZVJUST adjusts the history array on a change of order.
+C ZEWSET sets the error weight vector EWT before each step.
+C ZVNORM computes the weighted r.m.s. norm of a vector.
+C ZABSSQ computes the squared absolute value of a double complex z.
+C ZVSRCO is a user-callable routine to save and restore
+C the contents of the internal COMMON blocks.
+C ZACOPY is a routine to copy one two-dimensional array to another.
+C ZGEFA and ZGESL are routines from LINPACK for solving full
+C systems of linear algebraic equations.
+C ZGBFA and ZGBSL are routines from LINPACK for solving banded
+C linear systems.
+C DZSCAL scales a double complex array by a double prec. scalar.
+C DZAXPY adds a D.P. scalar times one complex vector to another.
+C ZCOPY is a basic linear algebra module from the BLAS.
+C DUMACH sets the unit roundoff of the machine.
+C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
+C error messages and warnings. XERRWD is machine-dependent.
+C Note: ZVNORM, ZABSSQ, DUMACH, IXSAV, and IUMACH are function routines.
+C All the others are subroutines.
+C The intrinsic functions called with double precision complex arguments
+C are: ABS, DREAL, and DIMAG. All of these are expected to return
+C double precision real values.
+C
+C-----------------------------------------------------------------------
+C
+C Type declarations for labeled COMMON block ZVOD01 --------------------
+C
+ DOUBLE PRECISION ACNRM, CCMXJ, CONP, CRATE, DRC, EL,
+ 1 ETA, ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU, TQ, TN, UROUND
+ INTEGER ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 1 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 2 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 3 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 4 NSLP, NYH
+C
+C Type declarations for labeled COMMON block ZVOD02 --------------------
+C
+ DOUBLE PRECISION HU
+ INTEGER NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+C Type declarations for local variables --------------------------------
+C
+ EXTERNAL ZVNLSD
+ LOGICAL IHIT
+ DOUBLE PRECISION ATOLI, BIG, EWTI, FOUR, H0, HMAX, HMX, HUN, ONE,
+ 1 PT2, RH, RTOLI, SIZE, TCRIT, TNEXT, TOLSF, TP, TWO, ZERO
+ INTEGER I, IER, IFLAG, IMXER, JCO, KGO, LENIW, LENJ, LENP, LENZW,
+ 1 LENRW, LENWM, LF0, MBAND, MFA, ML, MORD, MU, MXHNL0, MXSTP0,
+ 2 NITER, NSLAST
+ CHARACTER*80 MSG
+C
+C Type declaration for function subroutines called ---------------------
+C
+ DOUBLE PRECISION DUMACH, ZVNORM
+C
+ DIMENSION MORD(2)
+C-----------------------------------------------------------------------
+C The following Fortran-77 declaration is to cause the values of the
+C listed (local) variables to be saved between calls to ZVODE.
+C-----------------------------------------------------------------------
+ SAVE MORD, MXHNL0, MXSTP0
+ SAVE ZERO, ONE, TWO, FOUR, PT2, HUN
+C-----------------------------------------------------------------------
+C The following internal COMMON blocks contain variables which are
+C communicated between subroutines in the ZVODE package, or which are
+C to be saved between calls to ZVODE.
+C In each block, real variables precede integers.
+C The block /ZVOD01/ appears in subroutines ZVODE, ZVINDY, ZVSTEP,
+C ZVSET, ZVNLSD, ZVJAC, ZVSOL, ZVJUST and ZVSRCO.
+C The block /ZVOD02/ appears in subroutines ZVODE, ZVINDY, ZVSTEP,
+C ZVNLSD, ZVJAC, and ZVSRCO.
+C
+C The variables stored in the internal COMMON blocks are as follows:
+C
+C ACNRM = Weighted r.m.s. norm of accumulated correction vectors.
+C CCMXJ = Threshhold on DRC for updating the Jacobian. (See DRC.)
+C CONP = The saved value of TQ(5).
+C CRATE = Estimated corrector convergence rate constant.
+C DRC = Relative change in H*RL1 since last ZVJAC call.
+C EL = Real array of integration coefficients. See ZVSET.
+C ETA = Saved tentative ratio of new to old H.
+C ETAMAX = Saved maximum value of ETA to be allowed.
+C H = The step size.
+C HMIN = The minimum absolute value of the step size H to be used.
+C HMXI = Inverse of the maximum absolute value of H to be used.
+C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
+C HNEW = The step size to be attempted on the next step.
+C HRL1 = Saved value of H*RL1.
+C HSCAL = Stepsize in scaling of YH array.
+C PRL1 = The saved value of RL1.
+C RC = Ratio of current H*RL1 to value on last ZVJAC call.
+C RL1 = The reciprocal of the coefficient EL(1).
+C SRUR = Sqrt(UROUND), used in difference quotient algorithms.
+C TAU = Real vector of past NQ step sizes, length 13.
+C TQ = A real vector of length 5 in which ZVSET stores constants
+C used for the convergence test, the error test, and the
+C selection of H at a new order.
+C TN = The independent variable, updated on each step taken.
+C UROUND = The machine unit roundoff. The smallest positive real number
+C such that 1.0 + UROUND .ne. 1.0
+C ICF = Integer flag for convergence failure in ZVNLSD:
+C 0 means no failures.
+C 1 means convergence failure with out of date Jacobian
+C (recoverable error).
+C 2 means convergence failure with current Jacobian or
+C singular matrix (unrecoverable error).
+C INIT = Saved integer flag indicating whether initialization of the
+C problem has been done (INIT = 1) or not.
+C IPUP = Saved flag to signal updating of Newton matrix.
+C JCUR = Output flag from ZVJAC showing Jacobian status:
+C JCUR = 0 means J is not current.
+C JCUR = 1 means J is current.
+C JSTART = Integer flag used as input to ZVSTEP:
+C 0 means perform the first step.
+C 1 means take a new step continuing from the last.
+C -1 means take the next step with a new value of MAXORD,
+C HMIN, HMXI, N, METH, MITER, and/or matrix parameters.
+C On return, ZVSTEP sets JSTART = 1.
+C JSV = Integer flag for Jacobian saving, = sign(MF).
+C KFLAG = A completion code from ZVSTEP with the following meanings:
+C 0 the step was succesful.
+C -1 the requested error could not be achieved.
+C -2 corrector convergence could not be achieved.
+C -3, -4 fatal error in VNLS (can not occur here).
+C KUTH = Input flag to ZVSTEP showing whether H was reduced by the
+C driver. KUTH = 1 if H was reduced, = 0 otherwise.
+C L = Integer variable, NQ + 1, current order plus one.
+C LMAX = MAXORD + 1 (used for dimensioning).
+C LOCJS = A pointer to the saved Jacobian, whose storage starts at
+C WM(LOCJS), if JSV = 1.
+C LYH, LEWT, LACOR, LSAVF, LWM, LIWM = Saved integer pointers
+C to segments of ZWORK, RWORK, and IWORK.
+C MAXORD = The maximum order of integration method to be allowed.
+C METH/MITER = The method flags. See MF.
+C MSBJ = The maximum number of steps between J evaluations, = 50.
+C MXHNIL = Saved value of optional input MXHNIL.
+C MXSTEP = Saved value of optional input MXSTEP.
+C N = The number of first-order ODEs, = NEQ.
+C NEWH = Saved integer to flag change of H.
+C NEWQ = The method order to be used on the next step.
+C NHNIL = Saved counter for occurrences of T + H = T.
+C NQ = Integer variable, the current integration method order.
+C NQNYH = Saved value of NQ*NYH.
+C NQWAIT = A counter controlling the frequency of order changes.
+C An order change is about to be considered if NQWAIT = 1.
+C NSLJ = The number of steps taken as of the last Jacobian update.
+C NSLP = Saved value of NST as of last Newton matrix update.
+C NYH = Saved value of the initial value of NEQ.
+C HU = The step size in t last used.
+C NCFN = Number of nonlinear convergence failures so far.
+C NETF = The number of error test failures of the integrator so far.
+C NFE = The number of f evaluations for the problem so far.
+C NJE = The number of Jacobian evaluations so far.
+C NLU = The number of matrix LU decompositions so far.
+C NNI = Number of nonlinear iterations so far.
+C NQU = The method order last used.
+C NST = The number of steps taken for the problem so far.
+C-----------------------------------------------------------------------
+ COMMON /ZVOD01/ ACNRM, CCMXJ, CONP, CRATE, DRC, EL(13), ETA,
+ 1 ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU(13), TQ(5), TN, UROUND,
+ 3 ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 4 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 5 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 6 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 7 NSLP, NYH
+ COMMON /ZVOD02/ HU, NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+ DATA MORD(1) /12/, MORD(2) /5/, MXSTP0 /500/, MXHNL0 /10/
+ DATA ZERO /0.0D0/, ONE /1.0D0/, TWO /2.0D0/, FOUR /4.0D0/,
+ 1 PT2 /0.2D0/, HUN /100.0D0/
+C-----------------------------------------------------------------------
+C Block A.
+C This code block is executed on every call.
+C It tests ISTATE and ITASK for legality and branches appropriately.
+C If ISTATE .gt. 1 but the flag INIT shows that initialization has
+C not yet been done, an error return occurs.
+C If ISTATE = 1 and TOUT = T, return immediately.
+C-----------------------------------------------------------------------
+ IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601
+ IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602
+ IF (ISTATE .EQ. 1) GO TO 10
+ IF (INIT .NE. 1) GO TO 603
+ IF (ISTATE .EQ. 2) GO TO 200
+ GO TO 20
+ 10 INIT = 0
+ IF (TOUT .EQ. T) RETURN
+C-----------------------------------------------------------------------
+C Block B.
+C The next code block is executed for the initial call (ISTATE = 1),
+C or for a continuation call with parameter changes (ISTATE = 3).
+C It contains checking of all input and various initializations.
+C
+C First check legality of the non-optional input NEQ, ITOL, IOPT,
+C MF, ML, and MU.
+C-----------------------------------------------------------------------
+ 20 IF (NEQ .LE. 0) GO TO 604
+ IF (ISTATE .EQ. 1) GO TO 25
+ IF (NEQ .GT. N) GO TO 605
+ 25 N = NEQ
+ IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606
+ IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 607
+ JSV = SIGN(1,MF)
+ MFA = ABS(MF)
+ METH = MFA/10
+ MITER = MFA - 10*METH
+ IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608
+ IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608
+ IF (MITER .LE. 3) GO TO 30
+ ML = IWORK(1)
+ MU = IWORK(2)
+ IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609
+ IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610
+ 30 CONTINUE
+C Next process and check the optional input. ---------------------------
+ IF (IOPT .EQ. 1) GO TO 40
+ MAXORD = MORD(METH)
+ MXSTEP = MXSTP0
+ MXHNIL = MXHNL0
+ IF (ISTATE .EQ. 1) H0 = ZERO
+ HMXI = ZERO
+ HMIN = ZERO
+ GO TO 60
+ 40 MAXORD = IWORK(5)
+ IF (MAXORD .LT. 0) GO TO 611
+ IF (MAXORD .EQ. 0) MAXORD = 100
+ MAXORD = MIN(MAXORD,MORD(METH))
+ MXSTEP = IWORK(6)
+ IF (MXSTEP .LT. 0) GO TO 612
+ IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0
+ MXHNIL = IWORK(7)
+ IF (MXHNIL .LT. 0) GO TO 613
+ IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0
+ IF (ISTATE .NE. 1) GO TO 50
+ H0 = RWORK(5)
+ IF ((TOUT - T)*H0 .LT. ZERO) GO TO 614
+ 50 HMAX = RWORK(6)
+ IF (HMAX .LT. ZERO) GO TO 615
+ HMXI = ZERO
+ IF (HMAX .GT. ZERO) HMXI = ONE/HMAX
+ HMIN = RWORK(7)
+ IF (HMIN .LT. ZERO) GO TO 616
+C-----------------------------------------------------------------------
+C Set work array pointers and check lengths LZW, LRW, and LIW.
+C Pointers to segments of ZWORK, RWORK, and IWORK are named by prefixing
+C L to the name of the segment. E.g., segment YH starts at ZWORK(LYH).
+C Segments of ZWORK (in order) are denoted YH, WM, SAVF, ACOR.
+C Besides optional inputs/outputs, RWORK has only the segment EWT.
+C Within WM, LOCJS is the location of the saved Jacobian (JSV .gt. 0).
+C-----------------------------------------------------------------------
+ 60 LYH = 1
+ IF (ISTATE .EQ. 1) NYH = N
+ LWM = LYH + (MAXORD + 1)*NYH
+ JCO = MAX(0,JSV)
+ IF (MITER .EQ. 0) LENWM = 0
+ IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
+ LENWM = (1 + JCO)*N*N
+ LOCJS = N*N + 1
+ ENDIF
+ IF (MITER .EQ. 3) LENWM = N
+ IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
+ MBAND = ML + MU + 1
+ LENP = (MBAND + ML)*N
+ LENJ = MBAND*N
+ LENWM = LENP + JCO*LENJ
+ LOCJS = LENP + 1
+ ENDIF
+ LSAVF = LWM + LENWM
+ LACOR = LSAVF + N
+ LENZW = LACOR + N - 1
+ IWORK(17) = LENZW
+ LEWT = 21
+ LENRW = 20 + N
+ IWORK(18) = LENRW
+ LIWM = 1
+ LENIW = 30 + N
+ IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 30
+ IWORK(19) = LENIW
+ IF (LENZW .GT. LZW) GO TO 628
+ IF (LENRW .GT. LRW) GO TO 617
+ IF (LENIW .GT. LIW) GO TO 618
+C Check RTOL and ATOL for legality. ------------------------------------
+ RTOLI = RTOL(1)
+ ATOLI = ATOL(1)
+ DO 70 I = 1,N
+ IF (ITOL .GE. 3) RTOLI = RTOL(I)
+ IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
+ IF (RTOLI .LT. ZERO) GO TO 619
+ IF (ATOLI .LT. ZERO) GO TO 620
+ 70 CONTINUE
+ IF (ISTATE .EQ. 1) GO TO 100
+C If ISTATE = 3, set flag to signal parameter changes to ZVSTEP. -------
+ JSTART = -1
+ IF (NQ .LE. MAXORD) GO TO 200
+C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
+ CALL ZCOPY (N, ZWORK(LWM), 1, ZWORK(LSAVF), 1)
+ GO TO 200
+C-----------------------------------------------------------------------
+C Block C.
+C The next block is for the initial call only (ISTATE = 1).
+C It contains all remaining initializations, the initial call to F,
+C and the calculation of the initial step size.
+C The error weights in EWT are inverted after being loaded.
+C-----------------------------------------------------------------------
+ 100 UROUND = DUMACH()
+ TN = T
+ IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 110
+ TCRIT = RWORK(1)
+ IF ((TCRIT - TOUT)*(TOUT - T) .LT. ZERO) GO TO 625
+ IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO)
+ 1 H0 = TCRIT - T
+ 110 JSTART = 0
+ IF (MITER .GT. 0) SRUR = SQRT(UROUND)
+ CCMXJ = PT2
+ MSBJ = 50
+ NHNIL = 0
+ NST = 0
+ NJE = 0
+ NNI = 0
+ NCFN = 0
+ NETF = 0
+ NLU = 0
+ NSLJ = 0
+ NSLAST = 0
+ HU = ZERO
+ NQU = 0
+C Initial call to F. (LF0 points to YH(*,2).) -------------------------
+ LF0 = LYH + NYH
+ CALL F (N, T, Y, ZWORK(LF0), RPAR, IPAR)
+ NFE = 1
+C Load the initial value vector in YH. ---------------------------------
+ CALL ZCOPY (N, Y, 1, ZWORK(LYH), 1)
+C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
+ NQ = 1
+ H = ONE
+ CALL ZEWSET (N, ITOL, RTOL, ATOL, ZWORK(LYH), RWORK(LEWT))
+ DO 120 I = 1,N
+ IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621
+ 120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
+ IF (H0 .NE. ZERO) GO TO 180
+C Call ZVHIN to set initial step size H0 to be attempted. --------------
+ CALL ZVHIN (N, T, ZWORK(LYH), ZWORK(LF0), F, RPAR, IPAR, TOUT,
+ 1 UROUND, RWORK(LEWT), ITOL, ATOL, Y, ZWORK(LACOR), H0,
+ 2 NITER, IER)
+ NFE = NFE + NITER
+ IF (IER .NE. 0) GO TO 622
+C Adjust H0 if necessary to meet HMAX bound. ---------------------------
+ 180 RH = ABS(H0)*HMXI
+ IF (RH .GT. ONE) H0 = H0/RH
+C Load H with H0 and scale YH(*,2) by H0. ------------------------------
+ H = H0
+ CALL DZSCAL (N, H0, ZWORK(LF0), 1)
+ GO TO 270
+C-----------------------------------------------------------------------
+C Block D.
+C The next code block is for continuation calls only (ISTATE = 2 or 3)
+C and is to check stop conditions before taking a step.
+C-----------------------------------------------------------------------
+ 200 NSLAST = NST
+ KUTH = 0
+ GO TO (210, 250, 220, 230, 240), ITASK
+ 210 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
+ CALL ZVINDY (TOUT, 0, ZWORK(LYH), NYH, Y, IFLAG)
+ IF (IFLAG .NE. 0) GO TO 627
+ T = TOUT
+ GO TO 420
+ 220 TP = TN - HU*(ONE + HUN*UROUND)
+ IF ((TP - TOUT)*H .GT. ZERO) GO TO 623
+ IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
+ GO TO 400
+ 230 TCRIT = RWORK(1)
+ IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
+ IF ((TCRIT - TOUT)*H .LT. ZERO) GO TO 625
+ IF ((TN - TOUT)*H .LT. ZERO) GO TO 245
+ CALL ZVINDY (TOUT, 0, ZWORK(LYH), NYH, Y, IFLAG)
+ IF (IFLAG .NE. 0) GO TO 627
+ T = TOUT
+ GO TO 420
+ 240 TCRIT = RWORK(1)
+ IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
+ 245 HMX = ABS(TN) + ABS(H)
+ IHIT = ABS(TN - TCRIT) .LE. HUN*UROUND*HMX
+ IF (IHIT) GO TO 400
+ TNEXT = TN + HNEW*(ONE + FOUR*UROUND)
+ IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
+ H = (TCRIT - TN)*(ONE - FOUR*UROUND)
+ KUTH = 1
+C-----------------------------------------------------------------------
+C Block E.
+C The next block is normally executed for all calls and contains
+C the call to the one-step core integrator ZVSTEP.
+C
+C This is a looping point for the integration steps.
+C
+C First check for too many steps being taken, update EWT (if not at
+C start of problem), check for too much accuracy being requested, and
+C check for H below the roundoff level in T.
+C-----------------------------------------------------------------------
+ 250 CONTINUE
+ IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
+ CALL ZEWSET (N, ITOL, RTOL, ATOL, ZWORK(LYH), RWORK(LEWT))
+ DO 260 I = 1,N
+ IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510
+ 260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
+ 270 TOLSF = UROUND*ZVNORM (N, ZWORK(LYH), RWORK(LEWT))
+ IF (TOLSF .LE. ONE) GO TO 280
+ TOLSF = TOLSF*TWO
+ IF (NST .EQ. 0) GO TO 626
+ GO TO 520
+ 280 IF ((TN + H) .NE. TN) GO TO 290
+ NHNIL = NHNIL + 1
+ IF (NHNIL .GT. MXHNIL) GO TO 290
+ MSG = 'ZVODE-- Warning: internal T (=R1) and H (=R2) are'
+ CALL XERRWD (MSG, 50, 101, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG=' such that in the machine, T + H = T on the next step '
+ CALL XERRWD (MSG, 60, 101, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' (H = step size). solver will continue anyway'
+ CALL XERRWD (MSG, 50, 101, 1, 0, 0, 0, 2, TN, H)
+ IF (NHNIL .LT. MXHNIL) GO TO 290
+ MSG = 'ZVODE-- Above warning has been issued I1 times. '
+ CALL XERRWD (MSG, 50, 102, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' it will not be issued again for this problem'
+ CALL XERRWD (MSG, 50, 102, 1, 1, MXHNIL, 0, 0, ZERO, ZERO)
+ 290 CONTINUE
+C-----------------------------------------------------------------------
+C CALL ZVSTEP (Y, YH, NYH, YH, EWT, SAVF, VSAV, ACOR,
+C WM, IWM, F, JAC, F, ZVNLSD, RPAR, IPAR)
+C-----------------------------------------------------------------------
+ CALL ZVSTEP (Y, ZWORK(LYH), NYH, ZWORK(LYH), RWORK(LEWT),
+ 1 ZWORK(LSAVF), Y, ZWORK(LACOR), ZWORK(LWM), IWORK(LIWM),
+ 2 F, JAC, F, ZVNLSD, RPAR, IPAR)
+ KGO = 1 - KFLAG
+C Branch on KFLAG. Note: In this version, KFLAG can not be set to -3.
+C KFLAG .eq. 0, -1, -2
+ GO TO (300, 530, 540), KGO
+C-----------------------------------------------------------------------
+C Block F.
+C The following block handles the case of a successful return from the
+C core integrator (KFLAG = 0). Test for stop conditions.
+C-----------------------------------------------------------------------
+ 300 INIT = 1
+ KUTH = 0
+ GO TO (310, 400, 330, 340, 350), ITASK
+C ITASK = 1. If TOUT has been reached, interpolate. -------------------
+ 310 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
+ CALL ZVINDY (TOUT, 0, ZWORK(LYH), NYH, Y, IFLAG)
+ T = TOUT
+ GO TO 420
+C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
+ 330 IF ((TN - TOUT)*H .GE. ZERO) GO TO 400
+ GO TO 250
+C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
+ 340 IF ((TN - TOUT)*H .LT. ZERO) GO TO 345
+ CALL ZVINDY (TOUT, 0, ZWORK(LYH), NYH, Y, IFLAG)
+ T = TOUT
+ GO TO 420
+ 345 HMX = ABS(TN) + ABS(H)
+ IHIT = ABS(TN - TCRIT) .LE. HUN*UROUND*HMX
+ IF (IHIT) GO TO 400
+ TNEXT = TN + HNEW*(ONE + FOUR*UROUND)
+ IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
+ H = (TCRIT - TN)*(ONE - FOUR*UROUND)
+ KUTH = 1
+ GO TO 250
+C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
+ 350 HMX = ABS(TN) + ABS(H)
+ IHIT = ABS(TN - TCRIT) .LE. HUN*UROUND*HMX
+C-----------------------------------------------------------------------
+C Block G.
+C The following block handles all successful returns from ZVODE.
+C If ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
+C ISTATE is set to 2, and the optional output is loaded into the work
+C arrays before returning.
+C-----------------------------------------------------------------------
+ 400 CONTINUE
+ CALL ZCOPY (N, ZWORK(LYH), 1, Y, 1)
+ T = TN
+ IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
+ IF (IHIT) T = TCRIT
+ 420 ISTATE = 2
+ RWORK(11) = HU
+ RWORK(12) = HNEW
+ RWORK(13) = TN
+ IWORK(11) = NST
+ IWORK(12) = NFE
+ IWORK(13) = NJE
+ IWORK(14) = NQU
+ IWORK(15) = NEWQ
+ IWORK(20) = NLU
+ IWORK(21) = NNI
+ IWORK(22) = NCFN
+ IWORK(23) = NETF
+ RETURN
+C-----------------------------------------------------------------------
+C Block H.
+C The following block handles all unsuccessful returns other than
+C those for illegal input. First the error message routine is called.
+C if there was an error test or convergence test failure, IMXER is set.
+C Then Y is loaded from YH, and T is set to TN.
+C The optional output is loaded into the work arrays before returning.
+C-----------------------------------------------------------------------
+C The maximum number of steps was taken before reaching TOUT. ----------
+ 500 MSG = 'ZVODE-- At current T (=R1), MXSTEP (=I1) steps '
+ CALL XERRWD (MSG, 50, 201, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' taken on this call before reaching TOUT '
+ CALL XERRWD (MSG, 50, 201, 1, 1, MXSTEP, 0, 1, TN, ZERO)
+ ISTATE = -1
+ GO TO 580
+C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
+ 510 EWTI = RWORK(LEWT+I-1)
+ MSG = 'ZVODE-- At T (=R1), EWT(I1) has become R2 .le. 0.'
+ CALL XERRWD (MSG, 50, 202, 1, 1, I, 0, 2, TN, EWTI)
+ ISTATE = -6
+ GO TO 580
+C Too much accuracy requested for machine precision. -------------------
+ 520 MSG = 'ZVODE-- At T (=R1), too much accuracy requested '
+ CALL XERRWD (MSG, 50, 203, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' for precision of machine: see TOLSF (=R2) '
+ CALL XERRWD (MSG, 50, 203, 1, 0, 0, 0, 2, TN, TOLSF)
+ RWORK(14) = TOLSF
+ ISTATE = -2
+ GO TO 580
+C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
+ 530 MSG = 'ZVODE-- At T(=R1) and step size H(=R2), the error'
+ CALL XERRWD (MSG, 50, 204, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' test failed repeatedly or with abs(H) = HMIN'
+ CALL XERRWD (MSG, 50, 204, 1, 0, 0, 0, 2, TN, H)
+ ISTATE = -4
+ GO TO 560
+C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
+ 540 MSG = 'ZVODE-- At T (=R1) and step size H (=R2), the '
+ CALL XERRWD (MSG, 50, 205, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' corrector convergence failed repeatedly '
+ CALL XERRWD (MSG, 50, 205, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG = ' or with abs(H) = HMIN '
+ CALL XERRWD (MSG, 30, 205, 1, 0, 0, 0, 2, TN, H)
+ ISTATE = -5
+C Compute IMXER if relevant. -------------------------------------------
+ 560 BIG = ZERO
+ IMXER = 1
+ DO 570 I = 1,N
+ SIZE = ABS(ZWORK(I+LACOR-1))*RWORK(I+LEWT-1)
+ IF (BIG .GE. SIZE) GO TO 570
+ BIG = SIZE
+ IMXER = I
+ 570 CONTINUE
+ IWORK(16) = IMXER
+C Set Y vector, T, and optional output. --------------------------------
+ 580 CONTINUE
+ CALL ZCOPY (N, ZWORK(LYH), 1, Y, 1)
+ T = TN
+ RWORK(11) = HU
+ RWORK(12) = H
+ RWORK(13) = TN
+ IWORK(11) = NST
+ IWORK(12) = NFE
+ IWORK(13) = NJE
+ IWORK(14) = NQU
+ IWORK(15) = NQ
+ IWORK(20) = NLU
+ IWORK(21) = NNI
+ IWORK(22) = NCFN
+ IWORK(23) = NETF
+ RETURN
+C-----------------------------------------------------------------------
+C Block I.
+C The following block handles all error returns due to illegal input
+C (ISTATE = -3), as detected before calling the core integrator.
+C First the error message routine is called. If the illegal input
+C is a negative ISTATE, the run is aborted (apparent infinite loop).
+C-----------------------------------------------------------------------
+ 601 MSG = 'ZVODE-- ISTATE (=I1) illegal '
+ CALL XERRWD (MSG, 30, 1, 1, 1, ISTATE, 0, 0, ZERO, ZERO)
+ IF (ISTATE .LT. 0) GO TO 800
+ GO TO 700
+ 602 MSG = 'ZVODE-- ITASK (=I1) illegal '
+ CALL XERRWD (MSG, 30, 2, 1, 1, ITASK, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 603 MSG='ZVODE-- ISTATE (=I1) .gt. 1 but ZVODE not initialized '
+ CALL XERRWD (MSG, 60, 3, 1, 1, ISTATE, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 604 MSG = 'ZVODE-- NEQ (=I1) .lt. 1 '
+ CALL XERRWD (MSG, 30, 4, 1, 1, NEQ, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 605 MSG = 'ZVODE-- ISTATE = 3 and NEQ increased (I1 to I2) '
+ CALL XERRWD (MSG, 50, 5, 1, 2, N, NEQ, 0, ZERO, ZERO)
+ GO TO 700
+ 606 MSG = 'ZVODE-- ITOL (=I1) illegal '
+ CALL XERRWD (MSG, 30, 6, 1, 1, ITOL, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 607 MSG = 'ZVODE-- IOPT (=I1) illegal '
+ CALL XERRWD (MSG, 30, 7, 1, 1, IOPT, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 608 MSG = 'ZVODE-- MF (=I1) illegal '
+ CALL XERRWD (MSG, 30, 8, 1, 1, MF, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 609 MSG = 'ZVODE-- ML (=I1) illegal: .lt.0 or .ge.NEQ (=I2)'
+ CALL XERRWD (MSG, 50, 9, 1, 2, ML, NEQ, 0, ZERO, ZERO)
+ GO TO 700
+ 610 MSG = 'ZVODE-- MU (=I1) illegal: .lt.0 or .ge.NEQ (=I2)'
+ CALL XERRWD (MSG, 50, 10, 1, 2, MU, NEQ, 0, ZERO, ZERO)
+ GO TO 700
+ 611 MSG = 'ZVODE-- MAXORD (=I1) .lt. 0 '
+ CALL XERRWD (MSG, 30, 11, 1, 1, MAXORD, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 612 MSG = 'ZVODE-- MXSTEP (=I1) .lt. 0 '
+ CALL XERRWD (MSG, 30, 12, 1, 1, MXSTEP, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 613 MSG = 'ZVODE-- MXHNIL (=I1) .lt. 0 '
+ CALL XERRWD (MSG, 30, 13, 1, 1, MXHNIL, 0, 0, ZERO, ZERO)
+ GO TO 700
+ 614 MSG = 'ZVODE-- TOUT (=R1) behind T (=R2) '
+ CALL XERRWD (MSG, 40, 14, 1, 0, 0, 0, 2, TOUT, T)
+ MSG = ' integration direction is given by H0 (=R1) '
+ CALL XERRWD (MSG, 50, 14, 1, 0, 0, 0, 1, H0, ZERO)
+ GO TO 700
+ 615 MSG = 'ZVODE-- HMAX (=R1) .lt. 0.0 '
+ CALL XERRWD (MSG, 30, 15, 1, 0, 0, 0, 1, HMAX, ZERO)
+ GO TO 700
+ 616 MSG = 'ZVODE-- HMIN (=R1) .lt. 0.0 '
+ CALL XERRWD (MSG, 30, 16, 1, 0, 0, 0, 1, HMIN, ZERO)
+ GO TO 700
+ 617 CONTINUE
+ MSG='ZVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)'
+ CALL XERRWD (MSG, 60, 17, 1, 2, LENRW, LRW, 0, ZERO, ZERO)
+ GO TO 700
+ 618 CONTINUE
+ MSG='ZVODE-- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)'
+ CALL XERRWD (MSG, 60, 18, 1, 2, LENIW, LIW, 0, ZERO, ZERO)
+ GO TO 700
+ 619 MSG = 'ZVODE-- RTOL(I1) is R1 .lt. 0.0 '
+ CALL XERRWD (MSG, 40, 19, 1, 1, I, 0, 1, RTOLI, ZERO)
+ GO TO 700
+ 620 MSG = 'ZVODE-- ATOL(I1) is R1 .lt. 0.0 '
+ CALL XERRWD (MSG, 40, 20, 1, 1, I, 0, 1, ATOLI, ZERO)
+ GO TO 700
+ 621 EWTI = RWORK(LEWT+I-1)
+ MSG = 'ZVODE-- EWT(I1) is R1 .le. 0.0 '
+ CALL XERRWD (MSG, 40, 21, 1, 1, I, 0, 1, EWTI, ZERO)
+ GO TO 700
+ 622 CONTINUE
+ MSG='ZVODE-- TOUT (=R1) too close to T(=R2) to start integration'
+ CALL XERRWD (MSG, 60, 22, 1, 0, 0, 0, 2, TOUT, T)
+ GO TO 700
+ 623 CONTINUE
+ MSG='ZVODE-- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) '
+ CALL XERRWD (MSG, 60, 23, 1, 1, ITASK, 0, 2, TOUT, TP)
+ GO TO 700
+ 624 CONTINUE
+ MSG='ZVODE-- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) '
+ CALL XERRWD (MSG, 60, 24, 1, 0, 0, 0, 2, TCRIT, TN)
+ GO TO 700
+ 625 CONTINUE
+ MSG='ZVODE-- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) '
+ CALL XERRWD (MSG, 60, 25, 1, 0, 0, 0, 2, TCRIT, TOUT)
+ GO TO 700
+ 626 MSG = 'ZVODE-- At start of problem, too much accuracy '
+ CALL XERRWD (MSG, 50, 26, 1, 0, 0, 0, 0, ZERO, ZERO)
+ MSG=' requested for precision of machine: see TOLSF (=R1) '
+ CALL XERRWD (MSG, 60, 26, 1, 0, 0, 0, 1, TOLSF, ZERO)
+ RWORK(14) = TOLSF
+ GO TO 700
+ 627 MSG='ZVODE-- Trouble from ZVINDY. ITASK = I1, TOUT = R1. '
+ CALL XERRWD (MSG, 60, 27, 1, 1, ITASK, 0, 1, TOUT, ZERO)
+ GO TO 700
+ 628 CONTINUE
+ MSG='ZVODE-- ZWORK length needed, LENZW (=I1), exceeds LZW (=I2)'
+ CALL XERRWD (MSG, 60, 17, 1, 2, LENZW, LZW, 0, ZERO, ZERO)
+C
+ 700 CONTINUE
+ ISTATE = -3
+ RETURN
+C
+ 800 MSG = 'ZVODE-- Run aborted: apparent infinite loop '
+ CALL XERRWD (MSG, 50, 303, 2, 0, 0, 0, 0, ZERO, ZERO)
+ RETURN
+C----------------------- End of Subroutine ZVODE -----------------------
+ END
+*DECK ZVHIN
+ SUBROUTINE ZVHIN (N, T0, Y0, YDOT, F, RPAR, IPAR, TOUT, UROUND,
+ 1 EWT, ITOL, ATOL, Y, TEMP, H0, NITER, IER)
+ EXTERNAL F
+ DOUBLE COMPLEX Y0, YDOT, Y, TEMP
+ DOUBLE PRECISION T0, TOUT, UROUND, EWT, ATOL, H0
+ INTEGER N, IPAR, ITOL, NITER, IER
+ DIMENSION Y0(*), YDOT(*), EWT(*), ATOL(*), Y(*),
+ 1 TEMP(*), RPAR(*), IPAR(*)
+C-----------------------------------------------------------------------
+C Call sequence input -- N, T0, Y0, YDOT, F, RPAR, IPAR, TOUT, UROUND,
+C EWT, ITOL, ATOL, Y, TEMP
+C Call sequence output -- H0, NITER, IER
+C COMMON block variables accessed -- None
+C
+C Subroutines called by ZVHIN: F
+C Function routines called by ZVHIN: ZVNORM
+C-----------------------------------------------------------------------
+C This routine computes the step size, H0, to be attempted on the
+C first step, when the user has not supplied a value for this.
+C
+C First we check that TOUT - T0 differs significantly from zero. Then
+C an iteration is done to approximate the initial second derivative
+C and this is used to define h from w.r.m.s.norm(h**2 * yddot / 2) = 1.
+C A bias factor of 1/2 is applied to the resulting h.
+C The sign of H0 is inferred from the initial values of TOUT and T0.
+C
+C Communication with ZVHIN is done with the following variables:
+C
+C N = Size of ODE system, input.
+C T0 = Initial value of independent variable, input.
+C Y0 = Vector of initial conditions, input.
+C YDOT = Vector of initial first derivatives, input.
+C F = Name of subroutine for right-hand side f(t,y), input.
+C RPAR, IPAR = User's real/complex and integer work arrays.
+C TOUT = First output value of independent variable
+C UROUND = Machine unit roundoff
+C EWT, ITOL, ATOL = Error weights and tolerance parameters
+C as described in the driver routine, input.
+C Y, TEMP = Work arrays of length N.
+C H0 = Step size to be attempted, output.
+C NITER = Number of iterations (and of f evaluations) to compute H0,
+C output.
+C IER = The error flag, returned with the value
+C IER = 0 if no trouble occurred, or
+C IER = -1 if TOUT and T0 are considered too close to proceed.
+C-----------------------------------------------------------------------
+C
+C Type declarations for local variables --------------------------------
+C
+ DOUBLE PRECISION AFI, ATOLI, DELYI, H, HALF, HG, HLB, HNEW, HRAT,
+ 1 HUB, HUN, PT1, T1, TDIST, TROUND, TWO, YDDNRM
+ INTEGER I, ITER
+C
+C Type declaration for function subroutines called ---------------------
+C
+ DOUBLE PRECISION ZVNORM
+C-----------------------------------------------------------------------
+C The following Fortran-77 declaration is to cause the values of the
+C listed (local) variables to be saved between calls to this integrator.
+C-----------------------------------------------------------------------
+ SAVE HALF, HUN, PT1, TWO
+ DATA HALF /0.5D0/, HUN /100.0D0/, PT1 /0.1D0/, TWO /2.0D0/
+C
+ NITER = 0
+ TDIST = ABS(TOUT - T0)
+ TROUND = UROUND*MAX(ABS(T0),ABS(TOUT))
+ IF (TDIST .LT. TWO*TROUND) GO TO 100
+C
+C Set a lower bound on h based on the roundoff level in T0 and TOUT. ---
+ HLB = HUN*TROUND
+C Set an upper bound on h based on TOUT-T0 and the initial Y and YDOT. -
+ HUB = PT1*TDIST
+ ATOLI = ATOL(1)
+ DO 10 I = 1, N
+ IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
+ DELYI = PT1*ABS(Y0(I)) + ATOLI
+ AFI = ABS(YDOT(I))
+ IF (AFI*HUB .GT. DELYI) HUB = DELYI/AFI
+ 10 CONTINUE
+C
+C Set initial guess for h as geometric mean of upper and lower bounds. -
+ ITER = 0
+ HG = SQRT(HLB*HUB)
+C If the bounds have crossed, exit with the mean value. ----------------
+ IF (HUB .LT. HLB) THEN
+ H0 = HG
+ GO TO 90
+ ENDIF
+C
+C Looping point for iteration. -----------------------------------------
+ 50 CONTINUE
+C Estimate the second derivative as a difference quotient in f. --------
+ H = SIGN (HG, TOUT - T0)
+ T1 = T0 + H
+ DO 60 I = 1, N
+ 60 Y(I) = Y0(I) + H*YDOT(I)
+ CALL F (N, T1, Y, TEMP, RPAR, IPAR)
+ DO 70 I = 1, N
+ 70 TEMP(I) = (TEMP(I) - YDOT(I))/H
+ YDDNRM = ZVNORM (N, TEMP, EWT)
+C Get the corresponding new value of h. --------------------------------
+ IF (YDDNRM*HUB*HUB .GT. TWO) THEN
+ HNEW = SQRT(TWO/YDDNRM)
+ ELSE
+ HNEW = SQRT(HG*HUB)
+ ENDIF
+ ITER = ITER + 1
+C-----------------------------------------------------------------------
+C Test the stopping conditions.
+C Stop if the new and previous h values differ by a factor of .lt. 2.
+C Stop if four iterations have been done. Also, stop with previous h
+C if HNEW/HG .gt. 2 after first iteration, as this probably means that
+C the second derivative value is bad because of cancellation error.
+C-----------------------------------------------------------------------
+ IF (ITER .GE. 4) GO TO 80
+ HRAT = HNEW/HG
+ IF ( (HRAT .GT. HALF) .AND. (HRAT .LT. TWO) ) GO TO 80
+ IF ( (ITER .GE. 2) .AND. (HNEW .GT. TWO*HG) ) THEN
+ HNEW = HG
+ GO TO 80
+ ENDIF
+ HG = HNEW
+ GO TO 50
+C
+C Iteration done. Apply bounds, bias factor, and sign. Then exit. ----
+ 80 H0 = HNEW*HALF
+ IF (H0 .LT. HLB) H0 = HLB
+ IF (H0 .GT. HUB) H0 = HUB
+ 90 H0 = SIGN(H0, TOUT - T0)
+ NITER = ITER
+ IER = 0
+ RETURN
+C Error return for TOUT - T0 too small. --------------------------------
+ 100 IER = -1
+ RETURN
+C----------------------- End of Subroutine ZVHIN -----------------------
+ END
+*DECK ZVINDY
+ SUBROUTINE ZVINDY (T, K, YH, LDYH, DKY, IFLAG)
+ DOUBLE COMPLEX YH, DKY
+ DOUBLE PRECISION T
+ INTEGER K, LDYH, IFLAG
+ DIMENSION YH(LDYH,*), DKY(*)
+C-----------------------------------------------------------------------
+C Call sequence input -- T, K, YH, LDYH
+C Call sequence output -- DKY, IFLAG
+C COMMON block variables accessed:
+C /ZVOD01/ -- H, TN, UROUND, L, N, NQ
+C /ZVOD02/ -- HU
+C
+C Subroutines called by ZVINDY: DZSCAL, XERRWD
+C Function routines called by ZVINDY: None
+C-----------------------------------------------------------------------
+C ZVINDY computes interpolated values of the K-th derivative of the
+C dependent variable vector y, and stores it in DKY. This routine
+C is called within the package with K = 0 and T = TOUT, but may
+C also be called by the user for any K up to the current order.
+C (See detailed instructions in the usage documentation.)
+C-----------------------------------------------------------------------
+C The computed values in DKY are gotten by interpolation using the
+C Nordsieck history array YH. This array corresponds uniquely to a
+C vector-valued polynomial of degree NQCUR or less, and DKY is set
+C to the K-th derivative of this polynomial at T.
+C The formula for DKY is:
+C q
+C DKY(i) = sum c(j,K) * (T - TN)**(j-K) * H**(-j) * YH(i,j+1)
+C j=K
+C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, TN = TCUR, H = HCUR.
+C The quantities NQ = NQCUR, L = NQ+1, N, TN, and H are
+C communicated by COMMON. The above sum is done in reverse order.
+C IFLAG is returned negative if either K or T is out of bounds.
+C
+C Discussion above and comments in driver explain all variables.
+C-----------------------------------------------------------------------
+C
+C Type declarations for labeled COMMON block ZVOD01 --------------------
+C
+ DOUBLE PRECISION ACNRM, CCMXJ, CONP, CRATE, DRC, EL,
+ 1 ETA, ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU, TQ, TN, UROUND
+ INTEGER ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 1 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 2 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 3 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 4 NSLP, NYH
+C
+C Type declarations for labeled COMMON block ZVOD02 --------------------
+C
+ DOUBLE PRECISION HU
+ INTEGER NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+C Type declarations for local variables --------------------------------
+C
+ DOUBLE PRECISION C, HUN, R, S, TFUZZ, TN1, TP, ZERO
+ INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
+ CHARACTER*80 MSG
+C-----------------------------------------------------------------------
+C The following Fortran-77 declaration is to cause the values of the
+C listed (local) variables to be saved between calls to this integrator.
+C-----------------------------------------------------------------------
+ SAVE HUN, ZERO
+C
+ COMMON /ZVOD01/ ACNRM, CCMXJ, CONP, CRATE, DRC, EL(13), ETA,
+ 1 ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU(13), TQ(5), TN, UROUND,
+ 3 ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 4 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 5 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 6 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 7 NSLP, NYH
+ COMMON /ZVOD02/ HU, NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+ DATA HUN /100.0D0/, ZERO /0.0D0/
+C
+ IFLAG = 0
+ IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
+ TFUZZ = HUN*UROUND*SIGN(ABS(TN) + ABS(HU), HU)
+ TP = TN - HU - TFUZZ
+ TN1 = TN + TFUZZ
+ IF ((T-TP)*(T-TN1) .GT. ZERO) GO TO 90
+C
+ S = (T - TN)/H
+ IC = 1
+ IF (K .EQ. 0) GO TO 15
+ JJ1 = L - K
+ DO 10 JJ = JJ1, NQ
+ 10 IC = IC*JJ
+ 15 C = REAL(IC)
+ DO 20 I = 1, N
+ 20 DKY(I) = C*YH(I,L)
+ IF (K .EQ. NQ) GO TO 55
+ JB2 = NQ - K
+ DO 50 JB = 1, JB2
+ J = NQ - JB
+ JP1 = J + 1
+ IC = 1
+ IF (K .EQ. 0) GO TO 35
+ JJ1 = JP1 - K
+ DO 30 JJ = JJ1, J
+ 30 IC = IC*JJ
+ 35 C = REAL(IC)
+ DO 40 I = 1, N
+ 40 DKY(I) = C*YH(I,JP1) + S*DKY(I)
+ 50 CONTINUE
+ IF (K .EQ. 0) RETURN
+ 55 R = H**(-K)
+ CALL DZSCAL (N, R, DKY, 1)
+ RETURN
+C
+ 80 MSG = 'ZVINDY-- K (=I1) illegal '
+ CALL XERRWD (MSG, 30, 51, 1, 1, K, 0, 0, ZERO, ZERO)
+ IFLAG = -1
+ RETURN
+ 90 MSG = 'ZVINDY-- T (=R1) illegal '
+ CALL XERRWD (MSG, 30, 52, 1, 0, 0, 0, 1, T, ZERO)
+ MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) '
+ CALL XERRWD (MSG, 60, 52, 1, 0, 0, 0, 2, TP, TN)
+ IFLAG = -2
+ RETURN
+C----------------------- End of Subroutine ZVINDY ----------------------
+ END
+*DECK ZVSTEP
+ SUBROUTINE ZVSTEP (Y, YH, LDYH, YH1, EWT, SAVF, VSAV, ACOR,
+ 1 WM, IWM, F, JAC, PSOL, VNLS, RPAR, IPAR)
+ EXTERNAL F, JAC, PSOL, VNLS
+ DOUBLE COMPLEX Y, YH, YH1, SAVF, VSAV, ACOR, WM
+ DOUBLE PRECISION EWT
+ INTEGER LDYH, IWM, IPAR
+ DIMENSION Y(*), YH(LDYH,*), YH1(*), EWT(*), SAVF(*), VSAV(*),
+ 1 ACOR(*), WM(*), IWM(*), RPAR(*), IPAR(*)
+C-----------------------------------------------------------------------
+C Call sequence input -- Y, YH, LDYH, YH1, EWT, SAVF, VSAV,
+C ACOR, WM, IWM, F, JAC, PSOL, VNLS, RPAR, IPAR
+C Call sequence output -- YH, ACOR, WM, IWM
+C COMMON block variables accessed:
+C /ZVOD01/ ACNRM, EL(13), H, HMIN, HMXI, HNEW, HSCAL, RC, TAU(13),
+C TQ(5), TN, JCUR, JSTART, KFLAG, KUTH,
+C L, LMAX, MAXORD, N, NEWQ, NQ, NQWAIT
+C /ZVOD02/ HU, NCFN, NETF, NFE, NQU, NST
+C
+C Subroutines called by ZVSTEP: F, DZAXPY, ZCOPY, DZSCAL,
+C ZVJUST, VNLS, ZVSET
+C Function routines called by ZVSTEP: ZVNORM
+C-----------------------------------------------------------------------
+C ZVSTEP performs one step of the integration of an initial value
+C problem for a system of ordinary differential equations.
+C ZVSTEP calls subroutine VNLS for the solution of the nonlinear system
+C arising in the time step. Thus it is independent of the problem
+C Jacobian structure and the type of nonlinear system solution method.
+C ZVSTEP returns a completion flag KFLAG (in COMMON).
+C A return with KFLAG = -1 or -2 means either ABS(H) = HMIN or 10
+C consecutive failures occurred. On a return with KFLAG negative,
+C the values of TN and the YH array are as of the beginning of the last
+C step, and H is the last step size attempted.
+C
+C Communication with ZVSTEP is done with the following variables:
+C
+C Y = An array of length N used for the dependent variable vector.
+C YH = An LDYH by LMAX array containing the dependent variables
+C and their approximate scaled derivatives, where
+C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
+C j-th derivative of y(i), scaled by H**j/factorial(j)
+C (j = 0,1,...,NQ). On entry for the first step, the first
+C two columns of YH must be set from the initial values.
+C LDYH = A constant integer .ge. N, the first dimension of YH.
+C N is the number of ODEs in the system.
+C YH1 = A one-dimensional array occupying the same space as YH.
+C EWT = An array of length N containing multiplicative weights
+C for local error measurements. Local errors in y(i) are
+C compared to 1.0/EWT(i) in various error tests.
+C SAVF = An array of working storage, of length N.
+C also used for input of YH(*,MAXORD+2) when JSTART = -1
+C and MAXORD .lt. the current order NQ.
+C VSAV = A work array of length N passed to subroutine VNLS.
+C ACOR = A work array of length N, used for the accumulated
+C corrections. On a successful return, ACOR(i) contains
+C the estimated one-step local error in y(i).
+C WM,IWM = Complex and integer work arrays associated with matrix
+C operations in VNLS.
+C F = Dummy name for the user-supplied subroutine for f.
+C JAC = Dummy name for the user-supplied Jacobian subroutine.
+C PSOL = Dummy name for the subroutine passed to VNLS, for
+C possible use there.
+C VNLS = Dummy name for the nonlinear system solving subroutine,
+C whose real name is dependent on the method used.
+C RPAR, IPAR = User's real/complex and integer work arrays.
+C-----------------------------------------------------------------------
+C
+C Type declarations for labeled COMMON block ZVOD01 --------------------
+C
+ DOUBLE PRECISION ACNRM, CCMXJ, CONP, CRATE, DRC, EL,
+ 1 ETA, ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU, TQ, TN, UROUND
+ INTEGER ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 1 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 2 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 3 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 4 NSLP, NYH
+C
+C Type declarations for labeled COMMON block ZVOD02 --------------------
+C
+ DOUBLE PRECISION HU
+ INTEGER NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+C Type declarations for local variables --------------------------------
+C
+ DOUBLE PRECISION ADDON, BIAS1,BIAS2,BIAS3, CNQUOT, DDN, DSM, DUP,
+ 1 ETACF, ETAMIN, ETAMX1, ETAMX2, ETAMX3, ETAMXF,
+ 2 ETAQ, ETAQM1, ETAQP1, FLOTL, ONE, ONEPSM,
+ 3 R, THRESH, TOLD, ZERO
+ INTEGER I, I1, I2, IBACK, J, JB, KFC, KFH, MXNCF, NCF, NFLAG
+C
+C Type declaration for function subroutines called ---------------------
+C
+ DOUBLE PRECISION ZVNORM
+C-----------------------------------------------------------------------
+C The following Fortran-77 declaration is to cause the values of the
+C listed (local) variables to be saved between calls to this integrator.
+C-----------------------------------------------------------------------
+ SAVE ADDON, BIAS1, BIAS2, BIAS3,
+ 1 ETACF, ETAMIN, ETAMX1, ETAMX2, ETAMX3, ETAMXF, ETAQ, ETAQM1,
+ 2 KFC, KFH, MXNCF, ONEPSM, THRESH, ONE, ZERO
+C-----------------------------------------------------------------------
+ COMMON /ZVOD01/ ACNRM, CCMXJ, CONP, CRATE, DRC, EL(13), ETA,
+ 1 ETAMAX, H, HMIN, HMXI, HNEW, HRL1, HSCAL, PRL1,
+ 2 RC, RL1, SRUR, TAU(13), TQ(5), TN, UROUND,
+ 3 ICF, INIT, IPUP, JCUR, JSTART, JSV, KFLAG, KUTH,
+ 4 L, LMAX, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
+ 5 LOCJS, MAXORD, METH, MITER, MSBJ, MXHNIL, MXSTEP,
+ 6 N, NEWH, NEWQ, NHNIL, NQ, NQNYH, NQWAIT, NSLJ,
+ 7 NSLP, NYH
+ COMMON /ZVOD02/ HU, NCFN, NETF, NFE, NJE, NLU, NNI, NQU, NST
+C
+ DATA KFC/-3/, KFH/-7/, MXNCF/10/
+ DATA ADDON /1.0D-6/, BIAS1 /6.0D0/, BIAS2 /6.0D0/,
+ 1 BIAS3 /10.0D0/, ETACF /0.25D0/, ETAMIN /0.1D0/,
+ 2 ETAMXF /0.2D0/, ETAMX1 /1.0D4/, ETAMX2 /10.0D0/,
+ 3 ETAMX3 /10.0D0/, ONEPSM /1.00001D0/, THRESH /1.5D0/
+ DATA ONE/1.0D0/, ZERO/0.0D0/
+C
+ KFLAG = 0
+ TOLD = TN
+ NCF = 0
+ JCUR = 0
+ NFLAG = 0
+ IF (JSTART .GT. 0) GO TO 20
+ IF (JSTART .EQ. -1) GO TO 100
+C-----------------------------------------------------------------------
+C On the first call, the order is set to 1, and other variables are
+C initialized. ETAMAX is the maximum ratio by which H can be increased
+C in a single step. It is normally 10, but is larger during the
+C first step to compensate for the small initial H. If a failure
+C occurs (in corrector convergence or error test), ETAMAX is set to 1
+C for the next increase.
+C-----------------------------------------------------------------------
+ LMAX = MAXORD + 1
+ NQ = 1
+ L = 2
+ NQNYH = NQ*LDYH
+ TAU(1) = H
+ PRL1 = ONE
+ RC = ZERO
+ ETAMAX = ETAMX1
+ NQWAIT = 2
+ HSCAL = H
+ GO TO 200
+C-----------------------------------------------------------------------
+C Take preliminary actions on a normal continuation step (JSTART.GT.0).
+C If the driver changed H, then ETA must be reset and NEWH set to 1.
+C If a change of order was dictated on the previous step, then
+C it is done here and appropriate adjustments in the history are made.
+C On an order decrease, the history array is adjusted by ZVJUST.
+C On an order increase, the history array is augmented by a column.
+C On a change of step size H, the history array YH is rescaled.
+C-----------------------------------------------------------------------
+ 20 CONTINUE
+ IF (KUTH .EQ. 1) THEN
+ ETA = MIN(ETA,H/HSCAL)
+ NEWH = 1
+ ENDIF
+ 50 IF (NEWH .EQ. 0) GO TO 200
+ IF (NEWQ .EQ. NQ) GO TO 150
+ IF (NEWQ .LT. NQ) THEN
+ CALL ZVJUST (YH, LDYH, -1)
+ NQ = NEWQ
+ L = NQ + 1
+ NQWAIT = L
+ GO TO 150
+ ENDIF
+ IF (NEWQ .GT. NQ) THEN
+ CALL ZVJUST (YH, LDYH, 1)
+ NQ = NEWQ
+ L = NQ + 1
+ NQWAIT = L
+ GO TO 150
+ ENDIF
+C-----------------------------------------------------------------------
+C The following block handles preliminaries needed when JSTART = -1.
+C If N was reduced, zero out part of YH to avoid undefined references.
+C If MAXORD was reduced to a value less than the tentative order NEWQ,
+C then NQ is set to MAXORD, and a new H ratio ETA is chosen.
+C Otherwise, we take the same preliminary actions as for JSTART .gt. 0.
+C In any case, NQWAIT is reset to L = NQ + 1 to prevent further
+C changes in order for that many steps.
+C The new H ratio ETA is limited by the input H if KUTH = 1,
+C by HMIN if KUTH = 0, and by HMXI in any case.
+C Finally, the history array YH is rescaled.
+C-----------------------------------------------------------------------
+ 100 CONTINUE
+ LMAX = MAXORD + 1
+ IF (N .EQ. LDYH) GO TO 120
+ I1 = 1 + (NEWQ + 1)*LDYH
+ I2 = (MAXORD + 1)*LDYH
+ IF (I1 .GT. I2) GO TO 120
+ DO 110 I = I1, I2
+ 110 YH1(I) = ZERO
+ 120 IF (NEWQ .LE. MAXORD) GO TO 140
+ FLOTL = REAL(LMAX)
+ IF (MAXORD .LT. NQ-1) THEN
+ DDN = ZVNORM (N, SAVF, EWT)/TQ(1)
+ ETA = ONE/((BIAS1*DDN)**(ONE/FLOTL) + ADDON)
+ ENDIF
+ IF (MAXORD .EQ. NQ .AND. NEWQ .EQ. NQ+1) ETA = ETAQ
+ IF (MAXORD .EQ. NQ-1 .AND. NEWQ .EQ. NQ+1) THEN
+ ETA = ETAQM1
+ CALL ZVJUST (YH, LDYH, -1)
+ ENDIF
+ IF (MAXORD .EQ. NQ-1 .AND. NEWQ .EQ. NQ) THEN
+ DDN = ZVNORM (N, SAVF, EWT)/TQ(1)
+ ETA = ONE/((BIAS1*DDN)**(ONE/FLOTL) + ADDON)
+ CALL ZVJUST (YH, LDYH, -1)
+ ENDIF
+ ETA = MIN(ETA,ONE)
+ NQ = MAXORD
+ L = LMAX
+ 140 IF (KUTH .EQ. 1) ETA = MIN(ETA,ABS(H/HSCAL))
+ IF (KUTH .EQ. 0) ETA = MAX(ETA,HMIN/ABS(HSCAL))
+ ETA = ETA/MAX(ONE,ABS(HSCAL)*HMXI*ETA)
+ NEWH = 1
+ NQWAIT = L
+ IF (NEWQ .LE. MAXORD) GO TO 50
+C Rescale the history array for a change in H by a factor of ETA. ------
+ 150 R = ONE
+ DO 180 J = 2, L
+ R = R*ETA
+ CALL DZSCAL (N, R, YH(1,J), 1 )
+ 180 CONTINUE
+ H = HSCAL*ETA
+ HSCAL = H
+ RC = RC*ETA
+ NQNYH = NQ*LDYH
+C-----------------------------------------------------------------------
+C This section computes the predicted values by effectively
+C multiplying the YH array by the Pascal triangle matrix.
+C ZVSET is called to calculate all integration coefficients.
+C RC is the ratio of new to old values of the coefficient H/EL(2)=h/l1.
+C-----------------------------------------------------------------------
+ 200 TN = TN + H
+ I1 = NQNYH + 1
+ DO 220 JB = 1, NQ
+ I1 = I1 - LDYH
+ DO 210 I = I1, NQNYH
+ 210 YH1(I) = YH1(I) + YH1(I+LDYH)
+ 220 CONTINUE
+ CALL ZVSET
+ RL1 = ONE/EL(2)
+ RC = RC*(RL1/PRL1)
+ PRL1 = RL1
+C
+C Call the nonlinear system solver. ------------------------------------
+C
+ CALL VNLS (Y, YH, LDYH, VSAV, SAVF, EWT, ACOR, IWM, WM,
+ 1 F, JAC, PSOL, NFLAG, RPAR, IPAR)
+C
+ IF (NFLAG .EQ. 0) GO TO 450
+C-----------------------------------------------------------------------
+C The VNLS routine failed to achieve convergence (NFLAG .NE. 0).
+C The YH array is retracted to its values before prediction.
+C The step size H is reduced and the step is retried, if possible.
+C Otherwise, an error exit is taken.
+C-----------------------------------------------------------------------
+ NCF = NCF + 1
+ NCFN = NCFN + 1
+ ETAMAX = ONE
+ TN = TOLD
+ I1 = NQNYH + 1
+ DO 430 JB = 1, NQ
+ I1 = I1 - LDYH
+ DO 420 I = I1, NQNYH
+ 420 YH1(I) = YH1(I) - YH1(I+LDYH)
+ 430 CONTINUE
+ IF (NFLAG .LT. -1) GO TO 680
+ IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 670
+ IF (NCF .EQ. MXNCF) GO TO 670
+ ETA = ETACF
+ ETA = MAX(ETA,HMIN/ABS(H))
+ NFLAG = -1
+ GO TO 150
+C-----------------------------------------------------------------------
+C The corrector has converged (NFLAG = 0). The local error test is
+C made and control passes to statement 500 if it fails.
+C-----------------------------------------------------------------------
+ 450 CONTINUE
+ DSM = ACNRM/TQ(2)
+ IF (DSM .GT. ONE) GO TO 500
+C-----------------------------------------------------------------------
+C After a successful step, update the YH and TAU arrays and decrement
+C NQWAIT. If NQWAIT is then 1 and NQ .lt. MAXORD, then ACOR is saved
+C for use in a possible order increase on the next step.
+C If ETAMAX = 1 (a failure occurred this step), keep NQWAIT .ge. 2.
+C-----------------------------------------------------------------------
+ KFLAG = 0
+ NST = NST + 1
+ HU = H
+ NQU = NQ
+ DO 470 IBACK = 1, NQ
+ I = L - IBACK
+ 470 TAU(I+1) = TAU(I)
+ TAU(1) = H
+ DO 480 J = 1, L
+ CALL DZAXPY (N, EL(J), ACOR, 1, YH(1,J), 1 )
+ 480 CONTINUE
+ NQWAIT = NQWAIT - 1
+ IF ((L .EQ. LMAX) .OR. (NQWAIT .NE. 1)) GO TO 490
+ CALL ZCOPY (N, ACOR, 1, YH(1,LMAX), 1 )
+ CONP = TQ(5)
+ 490 IF (ETAMAX .NE. ONE) GO TO 560
+ IF (NQWAIT .LT. 2) NQWAIT = 2
+ NEWQ = NQ
+ NEWH = 0
+ ETA = ONE
+ HNEW = H
+ GO TO 690
+C-----------------------------------------------------------------------
+C The error test failed. KFLAG keeps track of multiple failures.
+C Restore TN and the YH array to their previous values, and prepare
+C to try the step again. Compute the optimum step size for the
+C same order. After repeated failures, H is forced to decrease
+C more rapidly.
+C-----------------------------------------------------------------------
+ 500 KFLAG = KFLAG - 1
+ NETF = NETF + 1
+ NFLAG = -2
+ TN = TOLD
+ I1 = NQNYH + 1
+ DO 520 JB = 1, NQ
+ I1 = I1 - LDYH
+ DO 510 I = I1, NQNYH
+ 510 YH1(I) = YH1(I) - YH1(I+LDYH)
+ 520 CONTINUE
+ IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 660
+ ETAMAX = ONE
+ IF (KFLAG .LE. KFC) GO TO 530
+C Compute ratio of new H to current H at the current order. ------------
+ FLOTL = REAL(L)
+ ETA = ONE/((BIAS2*DSM)**(ONE/FLOTL) + ADDON)
+ ETA = MAX(ETA,HMIN/ABS(H),ETAMIN)
+ IF ((KFLAG .LE. -2) .AND. (ETA .GT. ETAMXF)) ETA = ETAMXF
+ GO TO 150
+C-----------------------------------------------------------------------
+C Control reaches this section if 3 or more consecutive failures
+C have occurred. It is assumed that the elements of the YH array
+C have accumulated errors of the wrong order. The order is reduced
+C by one, if possible. Then H is reduced by a factor of 0.1 and
+C the step is retried. After a total of 7 consecutive failures,
+C an exit is taken with KFLAG = -1.
+C-----------------------------------------------------------------------
+ 530 IF (KFLAG .EQ. KFH) GO TO 660
+ IF (NQ .EQ. 1) GO TO 540
+ ETA = MAX(ETAMIN,HMIN/ABS(H))
+ CALL ZVJUST (YH, LDYH, -1)
+ L = NQ
+ NQ = NQ - 1
+ NQWAIT = L
+ GO TO 150
+ 540 ETA = MAX(ETAMIN,H