c
c latest revision date: august 30, 1993
c
c-----------------------------------------------------------------------
c    *** author ***
c
c author and contact
c
c                      luca dieci (*)
c                      school of mathematics
c                      georgia institute of technology
c                      atlanta, ga 30332-0160
c                      phone: (404) 853-9209
c                      e-mail: dieci@math.gatech.edu
c
c (*) please communicate promptly
c     any mistake you should find.
c
c-----------------------------------------------------------------------
c
      subroutine dresol(neq,neqlen,x,nx,y,ny,mf,t,tout,rtol,atol,
     *       itask,istate,iopt,rwork,lrex,iwork,liw,probl,rarr,iarr)
c
c Permission to use, copy, modify, and distribute this software for any
c purpose without fee is hereby granted, provided that this entire notice
c is included in all copies of any software which is or includes a copy
c or modification of this software and in all copies of the supporting
c documentation for such software.
c THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
c WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKE ANY
c REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
c OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
c
      external probl
      integer neq, neqlen, nx, ny, mf, itask, istate, iopt, lrex
      integer iwork, liw, iarr
      real x, y, t, tout, rtol, atol, rwork, rarr
      dimension neq(neqlen), iwork(liw), iarr(*)
      dimension x(nx,*), y(ny), rwork(lrex), rarr(*)
c----------------------------------------------------------------------
c     *** name ****
c we will refer to the full body of the programs here as the
c
c                   ----------------------
c                   |   DRESOL PACKAGE   |
c                   ----------------------
c
c-----------------------------------------------------------------------
c    *** purpose ***
c
c dresol solves the initial value problem for stiff or nonstiff
c systems of Matrix Differential Riccati Equations
c
c     dX/dt = A21(t) + A22(t)*X -X*A11(t) - X*A12(t)*X
c     + appropriate ICs.
c
c there are a number of options provided with this solver, for
c symmetric and unsymmetric riccati equations (see documentation
c below).  We will refer to the following partition of the underlying
c matrix A(t) in what follows.
c
c               nnq     nnp
c             |             |
c        nnq{ |A(t)  | A(t) |
c             | 11   |  12  |
c  A(t):=     | -----|----  |
c             |      |      |
c        nnp{ |A(t)  | A(t) |
c             | 21   |  22  |
c
c as usual, the purpose of the riccati equation is that one of
c decoupling the above matrix A(t).  In particular, to bring it
c to block upper triangular structure.  That is, if X(t) solves the
c riccati equation, then A(t) is transformed in the following
c
c             |~       ~    |
c             |A(t)  | A(t) |
c  ~          | 11   |  12  |          ~
c  A(t):=     | -----|----  |  , where A(t) = A(t)+A(t)*X(t)
c             |      | ~    |           11     11     12
c             |  0   | A(t) |      ~                      ~
c             |      |  22  |      A(t) = A(t)-X(t)*A(t), A(t) = A(t)
c                                   22     22        12    12     12
c
c in what follows, we will refer to the above partitionings.
c--------------------------------------------------------------------------
c    *** history ***
c
c the DRESOL package is a collection of subroutines for the
c numerical integration of matrix Riccati differential equations.
c the basic ODE solver on which this solver is based is the well-known
c stiff/nonstiff solver LSODE of Hindmarsh, written at
c LLL in the early '80s.  The modifications to the original
c code of Hindmarsh have been quite extensive, but
c nonetheless much of the original philosophy and heuristics for
c order selection, local error estimatation
c (hence stepsize change), and so forth, remain unchanged.
c
c
c subroutines dresol calls the subroutine drequs and they both
c constitute merely an interface between the user's driver and
c the differential riccati equation solver, subroutine drequ.
c This interface has been created to relieve the user of
c a tedious calling sequence and to keep the calling
c sequence to dresol conforming to the style of the lsode
c package of hindmarsh.
c
c
c-----------------------------------------------------------------------
c
c author: luca dieci, school of mathematics, Ga Tech
c         atlanta, ga 30332-0160
c         phone: (404) 853-9209
c         e-mail: dieci@math.gatech.edu
c
c     please address all correspondance related to the riccati
c     equation to this author, who is uniquely responsible
c     for the possible mulfunctioning of the
c     riccati solver.  please communicate promptly
c     any errors you should find.
cc
c-----------------------------------------------------------------------
c    *** references ***
c
c           Alan C. Hindmarsh, "Odepack, a systematized collection of
c           ode solvers", in Scientific Computing,
c           R. S. Stepleman et al. (eds.),
c           north-holland, Amsterdam, 1983, pp. 55-64.
c
c           Luca Dieci, "DRESOL. A differential riccati equations' solver:
c                        user's guide.  Version 1.0,"
c           Technical Report # 071090-52
c           School of Mathematics, Georgia Institute of Technology
c           Atlanta, GA 30332-0160
c
c           Luca Dieci, "Numerical integration of the differential
c           Riccati equation and some related issues,"
c           SIAM J. Numer. Analysis 29-3 (1992), pp. 781-815.
c
c------------------------------------------------------------------------
c    *** summary of usage ***
c
c communication between the user and the dresol package, for normal
c situations, is handled trough calls to the subroutines dresol.
c this routine perform some of the array splittings and arrays
c checkings. (other splittings and/or checkings performed in drequ).
c if the user desires, he might call drequ directly.
c we discourage this. In any case, here below we describe basically
c the full set of options available.  (see the description
c under the subroutine drequ for an extra avilable option on the type
c of error tolerances allowed, if needed).
c
c this summary indicates the typical sequence of operations required
c from the user in order to properly call dresol, hence drequs and drequ
c
c a. first provide a subroutine of the form..
c               subroutine probl (t,a,rarr,iarr)
c supplying the problem matrix in a
c
c
c b. next determine (or guess) whether or not the problem is stiff.
c stiffness occurs when the frechet derivative (jacobian) of the riccati
c equation has one of its eigenvalues with real part negative and
c large in magnitude, compared to the reciprocal of the t span of interest.
c if the problem is nonstiff, use a method flag  mf=10 or 20 (or mf=11,21)
c if it is stiff, we suggest the basic option provided by dresol
c which consists in using the backward differentiation formulas
c for carrying out one step of the integration, and perform a newton
c like iteration for the nolinear system.  when coupled with the chord
c method for the nonlinear iteration, this is the standard choice.
c this is obtained by requiring mf = 21 (and neq(5) = 0).  more choices
c are however available depending on the required strategy for solving the
c nonlinear system (the algebraic riccati equation) arising at each
c discretization step.  see the parameter neq(5) in the description below.
c
c
c c. write a main program which calls subroutine dresol 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 On the first call to drequ, supply arguments as explained in the
c discussion which follows.
c
c
c d. on return, the user needs to monitor the solution x (typically,
c computed at the requested tout point) and the status flag istate
c if the run has not been succesfull (as signalled by a negative value
c for istate, see below) the user might need to remove the causes for the
c insucces before calling dresol again.
c
c
c e. in case of a succesful run (positive istate), the user needs only to
c choose a new output point (tout) and call dresol again.
c
c
c N.B.: The user should set the value of the data ncpw in the routine
c       xerrwv accordingly to his own arithmetic environment,
c       and set the format of label 10 in subroutine xerrwv accordingly.
c       This ensures proper writing of error messages.
c       The default value at this moment is for computers with 4 bytes per
c       integer (ibm-360, ibm-370, ibm-303, ibm-43xx, ibm-pc, sun, sparc,..).
c --------------------------------------------------------------------------
c  *** interface ***
c
c the interface to dresol has been kept basically the same as
c the original interface of lsode for conformity and safety reasons.
c so doing, we have internally generated the necessary calling sequence
c for drequ.  the user of the dresol package should not attempt
c to call drequ directly, but call subroutine dresol.
c
c the interface to dresol consists of the following parts.
c
c i.   the call sequence to subroutine dresol, which is a driver
c      routine for the solver.  this includes descriptions of the call
c      sequence arguments and of internally supplied routines.
c      following these descriptions is a description of
c      optional inputs available through the call sequence, and then
c      a description of optional outputs (in the work arrays).
c
c ii.  descriptions of other routines in the dresol 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
c
c iii. descriptions of common blocks to be declared in overlay
c      or similar environments, or to be saved when doing an interrupt
c      of the problem and continued solution later.
c
c iv.  description of two routines in the dresol package, either of
c      which the user may replace with his own version, if desired.
c      these relate to the measurement of errors.
c      we recommend not to alter these routines unless absolutely
c      necessary as it might affect the performance of the method
c
c-----------------------------------------------------------------------
c **** part i.  call sequence. ***
c
c the call sequence parameters used for input only are
c     neq, neqlen, nx, ny, mf, tout, rtol, atol, itask, iopt,
c     lrex, liw, probl,
c and those used for both input and output are
c     x, (y), t, istate.
c the work arrays rwork and iwork are also used for conditional and
c optional inputs and optional outputs.  (the term output here refers
c to the return from subroutine dresol to the user-s calling program.)
c and the arrays rarr and iarr can also be used for communication
c by the user.
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 on input.
c
c -----------------------------------------------------------------------
c  *** explanation of parameters and basic options ***
c
c here below we provide a datailed explanation of all of the parameters in the
c calling sequence to dresol and of most of the available options.  a number
c of extra options are also available, but they should not be needed,
c and are summarized at the end of this block.
c
c probl  = name of subroutine giving the problem matrix a.
c          it must be declared external in the calling program,
c          and it must have the following form
c
c               subroutine probl (t,a,rarr,iarr)
c               real t, rarr, a
c               integer iarr
c               dimension a(nn,nn), rarr(*), iarr(*)
c	   where nn=nnp+nnq (see neq(2), neq(3))
c
c	   in this routine we must supply al of the entries of the
c          matrix a, regardless of its structure and of icontr (see neq(4))
c          we set a by loading a(i,j) with the (i,j)-th entry of
c          the (time dependent) matrix a.  the matrix a must be
c          be supplied entirely, even if some storage savings (due to
c          possible symmetries) had been possible.
c          The real and integer arrays rarr and iarr are provided
c          for possible communication between the user and the solver.
c          If the user decides to use them, then he needs to dimension
c          them oppurtunely in the calling program.  Otherwise, just
c          treat them as dummy arguments.
c         (N.B.: They could be possibly used for the reimbedding strategy.)
c
c neq    = integer array of dimension at least neqlen, used for
c          communication with the solver, whose entries are as follows:
c          neq(1)=neqn, this is defined as neqn=nnp*nnq, where
c          neq(2)=nnq, this is the dimension of the square matrix a11
c          neq(3)=nnp, this is the dimension of the square matrix a22
c          neq(4)=icontr, this must be set to 0 if the riccati equation is
c                 of the type arising in optimal control from a hamiltonian
c                 matrix with guaranteed existence and stability properties.
c                 in this case, explicit exploitation of the structure
c                 is performed
c          neq(5)=nonlin, this parameter specifies how the nonlinear
c                 system arising at each step should be solved in case the
c                 user has chosen mf=11 or 21.  allowed values are
c                 nonlin=0,1 where 0 corresponds to the chord method, and
c                 1 corresponds to the quasi-newton method. a value of 1 is
c                 recommended when more accurate eigenvalue information
c                 is desired to possibly monitor stability or the like.
c                 usually, nonlin=0 is a cheaper choice (the default).
c          neq(6)=nsuper, this parameter specifies whether we want (nsuper=1)
c                 or not (nsuper=0) to perform the superstability check.
c                 it makes sense only for mf=11,21, because the test relies
c                 on eigenvalues of the jacobian, and we further suggest
c                 to use it only with mf=21.
c                 This test tries to avoid that with mf=21 (or 11)
c                 the solver steps on a wrong solution branch, missing a
c                 transition region.  The user must posses a complete
c                 understanding of the implications of using nsuper=1, if
c                 he is going to use this value (see the ref.s).  Otherwise,
c                 we recommend to let nsuper=0 (the default) and
c                 monitor iwork(2) (see optional output).
c                 With nsuper=1, the code attempts to avoid
c                 superstability occurrences both for neq(5)=0,1 but
c                 it works much better for neq(5)=1, because the code has
c                 current eigenvalues information.  We suggest
c                 to use nsuper=1 only with neq(5)=1.
c                 The superstability test performs satisfactorily
c                 when (i) there is a stable trajectory for the solution
c                 and the solution is (initially) along this trajectory
c                 (ii) the solution is along an unstable trajectory and we
c                 want to be cautious.  Again, we refer to the papers given
c                 in the introduction for a more detailed explanation.
c           N.B.: the solver tries to detect/avoid the superstability
c                 occurrences (see the documentation in the subroutine stode)
c                 on a given step, up to 6 step reductions, and unless the
c                 stepsize has become so small to lead to a loss of meaning.
c                 if the code fails to recover (due to these occurrences)
c                 from the possible superstability, then it does not try
c                 any more to recover for the given run and it signals its
c                 failing by setting nsuper=0, so that the user should
c                 monitor this quantity if he originally set nsuper=1.
c                 (the code prints a message saying that this occurred).
c
c neqlen = dimension of neq.  it must have neqlen.ge.6.
c
c x      = riccati matrix of dimension (nx,nnq) containing the initial
c          values of the riccati solution.  user must supply all entries of
c          x even if icontr=0, and hence there could have been symmetries in x
c          (which will be exploited internally by the code, in this case)
c          on output, x contains the matrix of computed values x(t).
c
c nx     = row dimension of the matrix x, must have nx.ge.nnp
c
c y      = work array of length ny used internally by the code, whose
c          content is usually of no relevance to the user.  the first neqn
c          components of y contain the riccati matrix x unrolled
c          according to y(k)=x(i,j), k=nnp*(j-1)+i, i=1,nnp, j=1,nnq.
c          the array y is passed as the y argument in all calls.
c          y(neqn+1),... are used to store the problem matrix
c          and its partitioned form as needed.
c
c ny	 = dimension of y of length at least neqn+2*nn*nn
c
c t      = on input, the initial value of the independent variable.
c          on output, value of independent variable (normally tout),
c          where the solution has been evaluated.  on an error return,
c          t is simply the farthest point reached.
c
c tout   = first point where output is desired, used only for input.
c          notice that we allow for either t.gt.tout or t.lt.tout
c          and so both forward and backward integration are allowed.
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.
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          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 outputs, below, for
c          tcur and hu).
c
c rtol   = relative tolerance parameter (scalar).
c
c atol   = absolute tolerance parameter (scalar).
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
c          thus the local error test passes if, in each component,
c          either the absolute error is less than atol,
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 for pure relative error control (not to use if
c          some solution component is known to vanish).
c          the values of rtol and atol should all be non-negative.
c          if the above choice (with atol, rtol scalar) is not suitable,
c          more general choices can be obtained.  we refer to the
c          supplementary documentation at the beginning of the
c          subroutine drequ below.
c          caution.. actual (global) errors may exceed these
c          local tolerances, so choose them conservatively.
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 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 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          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          state of the calculation.
c          on 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(4), neq(5), neq(6),
c             rtol, atol, iopt, lrex, liw, mf,
c             and any of the optional inputs except h0.
c             (of course, above changes can be made subject to mutual
c             compatibilities).
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 for the
c          purpose of outputting the initial conditions.)
c          thus the first call for which tout .ne. t requires
c          istate = 1 on input.
c          on output, istate has the following values and meanings.
c           1  means nothing was done, as tout was equal to t with
c              istate = 1 on input.  (however, an internal counter was
c              set to detect and prevent repeated calls of this type.)
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 below on optional inputs).
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=0.0)
c              was requested, but a variable has now vanished.
c              the integration was successful as far as t.
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 inputs, before
c          calling the solver again.
c
c iopt   = an integer flag to specify whether or not any optional
c          inputs are being used on this call.  input only.
c          the optional inputs are listed separately below.
c          iopt = 0 (the default) means no optional inputs are being used.
c                   default values will be used in all cases.
c          iopt = 1 means one or more optional inputs are being used.
c
c rwork  = a real work array of length at least lrex
c          the content of rwork can be very relevant to the user,
c          and some of it is explained in the optional output later on.
c          the first 20 words of rwork are reserved for conditional
c          and optional inputs and optional outputs.
c
c lrex   = declared length of rwork; we must have lrex at least
c             20 + neqn*(maxord + 4) + lns    where
c          neqn = neq(1)
c          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
c                   smaller value is given as an optional input),
c          lns = 0                           if miter = 0,
c          lns = 2+2*(nnp**2+nnq**2)+3*nn    if miter = 1
c          regardless of nonlin=0,1.
c          (see the mf description for meth and miter.)
c          this paramter lrex will be checked by the solver.
c          thus if maxord has its default value this length is..
c             20 + 16*neqn                  for mf = 10,
c             20 + 16*neqn + lns            for mf = 11
c             20 +  9*neqn                  for mf = 20,
c             20 +  9*neqn + lns            for mf = 21
c note..  the work arrays must not be altered between calls to drequ
c for the same problem.
c
c iwork  = integer work array of length at least liw
c          the first few words of iwork are used for conditional
c          and optional input and optional output
c
c liw    = declared length of iwork; we must have liw at least
c             20        for mf = 10, 20
c             20 + nn   for mf = 11, 21
c          this parameter liw will be checked by the solver
c
c mf     = the method flag.  used only for input.  the legal values of
c          mf are 10, 11, 20, 21.
c          mf has decimal digits meth and miter.. mf = 10*meth + miter.
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 an internally generated
c                      jacobian, kept in compact matrix form.
c          if miter = 1, the user must also specify the particular
c          choice of nonlinear solver he desires, as explained under
c          the parameter neq(5).
c          To summarize, the following are available and suggested options
c          10 for nonstiff (adams) method, no jacobian used, simple
c             functional fixed point iteration for nonlinear system
c          11 for adams method, but with jacobian used in a newton type
c             iteration (see neq(5))
c          20 for nonstiff (bdf) method, functional fixed point iteration
c          21 for stiff (bdf) method, but with jacobian used (recommended
c             value for bdf method).
c          We recommend to use mf=10 or 20 for nonstiff Riccati
c          equations and mf=21 for stiif equations.
c
c rarr,  = real and integer arrays which can be used by the user for
c   iarr     communication between the calling program and routines
c          probl (hence all of those called by probl).  dresol does not
c          touch these arrays, but it passes them.  If the user does
c          not need them ,just treat them as dummy arguments.  if user
c          decides to use them, they must be declared as arrays of
c          appropriate length in the calling program and in probl.
c          they can be put to use if user wants to perform a reimbedding
c          strategy, in the usual ways for riccati equations.
c
c note that the main program must declare arrays neq, x, y, rwork, iwork, and
c rarr and iarr if they are used.
c
c-----------------------------------------------------------------------
c  *** optional inputs ***
c
c the following is a list of the optional inputs 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 these inputs requires iopt = 1, and in that
c case all of these inputs are examined.  a value of zero for any
c of these optional inputs will cause the default value to be used.
c thus to use a subset of the optional inputs, 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 tcrit   rwork(1)  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 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               note: the default value is 0.0, so be careful.
c                   (this lower bound is not enforced on the final
c                   step before reaching tcrit 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 *** optional outputs ***
c
c as optional additional output from dresol, the variables listed
c below are quantities related to the performance of dresol
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 these outputs are defined
c on any successful return from drequ, 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, lenrw, and leniw.
c on any error return, outputs 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.  on output, tcur
c                   will always be at least as far as the argument
c                   t, 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).
c                   if 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 iwm(1)   iwork(1)  if miter=1 (mf=11 or 21), this entry is set to
c                    1 if during integration the solver found purely
c                    imaginary eigenvalues.  it is set to 0 otherwise.
c                    of course, this is not an error condition, but it
c                    might indicate that the user might want to use
c                    the adams formulas if he is using the backward-
c                    differentiation formulas.
c
c iwm(2)   iwork(2)  if miter=1 (mf=11 or 21), this entry os set to
c                    1 if during integration the solver found unstable
c                    eigenvalues with respect to the direction of
c                    integration.  it is set to 0 otherwise.
c                    this is not an error condition, but the user should
c                    monitor the eigenvalues to be certain that the
c                    integrator is performing adequately.
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                   in this code, every jacobian evaluation corresponds
c                   to the solution of a sylvester equation, hence
c                   to a factorization of the blocks
c                      ~    ~
c                      A  , A   via the bartels-stewart algorithm
c                       11   22
c                   (see the documentation for prepj).
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 lenrw   iwork(17) 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(18) 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
c the following arrays are segments of the rwork array which
c can be of interest to the user as optional outputs.
c for each array, the table below gives its internal name,
c its base address in rwork, and its description.
c
c
c N.B.: Important eigenvalue information is in these arrays in the entries
c       era, erb, eia, aib, itypea, itypeb, as described below.
c
c name    base address      description
c
c in the rwork array
c
c yh      21             the nordsieck history array, of size neqn by
c                        (nqcur + 1). 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 solution,
c                        evaluated at t = tcur.
c
c era      lwm+nttp      only if miter=1 (mf=11,21), and we have set
c                        lwm=20+(maxord+1)*neqn, nttp=3+2*nnp**2+2*nnq**2+nnp.
c                        array of dimension nnp, of real parts of ev.s
c                    --> if icontr=0, it contains the real parts of the
c                        eigenvalues of the closed loop spectrum matrix
c                        ~  ^            ^
c                        A (t   ), where t    is the farthest point
c                         11 n+1          n+1
c                        reached in the integration (not necessarily
c                        the same as tout), when nonlin=1, and otherwise
c                        is the last point where the code performed a
c                        jacobian evaluation/factorization if nonlin=0.
c                        These eigenvalues are ordered from smallest
c                        to largest real parts
c                    --> if icontr=1, it contains the real parts of the
c                        eigenvalues of the (nnp*nnp) lower matrix
c                        ~  ^            ^
c                        A (t   ), where t    is as above.
c                         22 n+1          n+1
c                        The eigenvalues are ordered from smallest to
c                        largest real part.
c
c eia      lwm+ntt2p     only when miter=1, and we set ntt2p=nttp+nnp.
c                        array of dimension nnp, containing the imaginary
c                        parts of the eigenvalues whose real parta are era
c                        (see aabove entry), and they are stored as in the
c                        way explained in era, that is eia(i) is imaginary
c                        part corresponding to era(i).
c
c erb      lwm+ntt1q     only when miter=1 and icontr=1, and we have
c                        set ntt1q=lwm+nttp+2*nnp+nnq.
c                        array of dimension nnq, containing the real
c                        parts of the eigenvalues of the upper matrix
c                        ~  ^            ^
c                        A (t   ), where t    is as above.
c                         11 n+1          n+1
c                        The eigenvalues are ordered from smallest to
c                        largest real part.
c
c eib       lwm+ntt2q    only when miter=1 and icontr=1, where
c                        we set ntt2q=ntt1q+nnq.
c                        array of dimension nnq, containing the imaginary
c                        parts eib(i), of the eigenvalues whose real parts
c                        are given by erb(i).
c
c acor     lenrw-neq+1   array of size neqn used for the accumulated
c                        corrections on each step, scaled on output
c                        to represent the estimated local error in x (i.e.
c                        in y) 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 dresol.
c
c in the iwork array
c
c itypea   21            array of size nnp, whose entries specify whether
c                        the i-th eigenvalue corresponding to (era(i),eia(i))
c                        is real, purely imaginary, or complex with nonzero
c                        real part.  see the documnetation to the
c                        subroutines axpxb, atxpxa, hqr3.
c
c itypeb   21+nnp        array of size nnq, whose entries are as for
c                        itypea, but now they are relative to the
c                        eigenvalues (erb(i),eib(i)).
c
c N.B.: further information which might be of some use is
c       given by the orthogonal matrices reducing the sylvester equation
c       blocks, when miter=1.  if interested, see description in prepj
c       and connected routines.
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 dresol.
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 dresol, 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 dresol.
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 srcom(rsav,isav,job) saves and restores the contents of
c                             the internal common blocks used by
c                             dresol (see part iii below).
c                             rsav must be a real array of length 222
c                             or more, and isav must be an integer
c                             array of length 66 or more.
c                             job=1 means save common into rsav/isav.
c                             job=2 means restore common from rsav/isav.
c                                srcom is useful if one is
c                             interrupting a run and restarting
c                             later, or alternating between two or
c                             more problems solved with dresol.
c
c   call intdy(,,,,,)         provide derivatives of x (i.e. of y), of
c        (see below)          various orders, at a specified point t, if
c                             desired.  it may be called only after
c                             a successful return from dresol.
c
c the detailed instructions for using intdy are as follows.
c the form of the call is..
c
c   call intdy (t, k, rwork(21), neqn, 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 drequ).
c             for valid results, t must lie between tcur - hu and tcur.
c             (see optional outputs 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 outputs).  the capability corresponding
c             to k = 0, i.e. computing y(t), is already provided
c             by drequ directly.  since nqcur .ge. 1, the first
c             derivative dx/dt is always available with intdy.
c rwork(21) = the base address of the history array yh.
c neqn      = column length of yh, equal to the value of neq(1).
c
c the output parameters are..
c
c dky       = a real array of length neqn containing the computed value
c             of the k-th derivative of the riccati solution x(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
c if dresol 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 dresol,
c   (2) the three internal common blocks
c         /ls0001/  of length  257  (218 single precision words
c                         followed by 39 integer words),
c         /eh0001/  of length  2 (integer words),
c         /intdre/ of length 29 (4 single precision words
c                         followed by 25 integer words).
c
c if dresol 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 three common blocks in his main program to insure
c that their contents are preserved.
c
c if the solution of a given problem by dresol is to be interrupted
c and then later continued, such as when restarting an interrupted run
c or alternating between two or more problems, the user should save,
c following the return from the last dresol call prior to the
c interruption, the contents of the call sequence variables and the
c internal common blocks, and later restore these values before the
c next dresol call for that problem.  to save and restore the common
c blocks, use subroutine srcom (see part ii above).
c
c-----------------------------------------------------------------------
c *** part iv.  optionally replaceable solver routines ***
c
c below are descriptions of two routines in the dresol 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) ewset.
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 rtol atol above.. (see documentation of drequ if
c need further atol/rtol possibilities, such as vector tolerances)
c     subroutine ewset (neqn, itol, rtol, atol, ycur, ewt)
c where neqn, rtol, and atol are as in the dresol call sequence,
c itol = 1 in this version (unless user decides to call drequ directly,
c and set itol to a different value, but this is discouraged),
c ycur contains the current dependent variable vector, and
c ewt is the array of weights set by ewset.
c
c if the user supplies this subroutine, it must return in ewt(i)
c (i = 1,...,neqn) a positive quantity suitable for comparing errors
c in y(i) to (recall that the riccati matrix x is stored internally
c in the vector y, columnwise).
c the ewt array returned by ewset is passed to the vnorm
c routine (see below), and also used by dresol in the computation
c of the optional output imxer.
c
c in the user-supplied version of ewset, it may be desirable to use
c the current values of derivatives of the riccati solution x.
c Derivatives up to order nq are available from the history array yh,
c described above under optional outputs.
c In ewset, yh is identical to the ycur array,
c extended to nq + 1 columns with a column length of neqn 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 the quantities nq, neqn, h, and nst can be obtained by including
c in ewset the statements..
c     common /ls0001/ rls(218),ils(39)
c     nq = ils(35)
c     neqn = ils(14)
c     nst = ils(36)
c     h = rls(212)
c thus, for example, the current value of dy/dt can be obtained as
c ycur(neqn+i)/h  (i=1,...,neq)  (and the division by h is
c unnecessary when nst = 0).
c
c (b) vnorm.
c the following is a real function routine which computes the weighted
c max-norm norm of a vector v..
c     d = vnorm (n, v, w)
c where..
c   n = the length of the vector,
c   v = real array of length n containing the vector,
c   w = real array of length n containing weights,
c   d = max ( abs(v(i)*w(i)), i=1, n )
c vnorm is called with n = neqn and with w(i) = 1.0/ewt(i), where
c ewt is as set by subroutine ewset.
c
c if the user supplies this function, it should return a non-negative
c value of vnorm suitable for use in the error control in drequ.
c none of the arguments should be altered by vnorm.
c for example, a user-supplied vnorm routine might..
c   -substitute the root-mean-square norm for the max-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 *** other routines in the dresol package ***
c
c in addition to subroutine dresol, drequs and drequ, which effectively
c perform initializations, decide upon returns-errors, and so forth,
c in such a way to keep the code as close to lsode as possible,
c the dresol package also includes the following subroutines and
c function routines.. (some of these are straight out of the lsode
c package, other are modified from lsode, and other are original to dresol).
c  intdy    computes an interpolated value of the y vector (which is
c           the unrolled riccati matrix x) at t = tout.
c  stode    is the core integrator, which does one step of the
c           integration and the associated error control
c           and eigenvalues checking.  if nsuper=1 also the superstability
c           testing is performed here.
c  cfode    sets all method coefficients and test constants.
c  f, doblks, and foftex compute the right hand side of the riccati
c           equation by maximizing efficiency (also depending on icontr)
c  prepj, jac  manage computation and preprocess of the jacobian matrix and
c           the newton iteration matrix p = i - h*l0*j.  there are several
c           important technical details here.  these routines call
c           the routines axpxbd, atxxad, orthes, ortran, hqr3, exchng,
c           qrstep, split, to perform an ordered eigenvalue decomposition.
c  solsy    manages solution of linear system in chord iteration and in the
c           quasi-newton iteration (nonlin=0,1).  it calls the routines
c           axpxbs, shrslv, atxxas, symslv, sysslv, to solve the decomposed
c           sylvester equations.
c  ewset    sets the error weight vector ewt before each step.
c  vnorm    computes the weighted max norm of a vector.
c  srcom    is a user-callable routine to save and restore
c           the contents of the internal common blocks.
c  r1mach   computes the unit roundoff in a machine-independent manner.
c  xerrwv, xsetun, and xsetf   handle the printing of all error
c           messages and warnings.  xerrwv is machine-dependent.
c note..  vnorm and r1mach are function routines.
c all the others are subroutines.
c
c there are also two block data in dresol, atlan1 and atlan2 which set
c the constant illin, ntrep (in atlan1) and msflg, lunit (in atlan2)
c (see part ii. other routines callable, for replacing these entries).
c
c the intrinsic and external routines used by dresol are..
c abs, amax1, amin1, float, max0, min0, mod, sign, sqrt, and write.
c *** end of prologue ***
      integer neqn, ney
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      real reuvct, reuvpt, tcj, tpj
c-----------------------------------------------------------------------
c the following internal common block contains
c (a) variables which are communicated between subroutines.  these are
c variables needed for arraies spitting and other checkings.
c the structure of the block is as follows..  all real variables are
c listed first, followed by all integers.  the block is declared in
c subroutines dresol, drequs, drequ, stode, f, prepj, and solsy.
c most of the variables of this block are of no interest to the
c user, but he might find interesting the quantity tcj and tpj:
c they are the last (and the previous to the last) time at which
c a jacobian was decomposed/updated.
c-----------------------------------------------------------------------
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
c here we perform an array splitting and check for the first time
c some of the arrays dimensions
c (further checking will be done in drequ)
c the checking is done only for first call or if input has been changed
      if (istate.ne.1.and.istate.ne.3) goto 10
      neqn=neq(1)
      nnq=neq(2)
      nnp=neq(3)
      icontr=neq(4)
      if (icontr .ne. 0) icontr = 1
      if (icontr.eq.0.and.nnq.ne.nnp) then
         istate=-3
         call xerrwv(45hdresol--  wrong icontr-nnp-nnq: run aborted  ,
     1      45,1,2,0,0,0,0,0.e0,0.e0)
         return
      endif
      nonlin=neq(5)
      if (nonlin.ne.1) nonlin=0
      nsuper=neq(6)
      if (nsuper.ne.1) nsuper=0
      nn=nnp+nnq
      if (neqlen.lt.6) then
         istate=-3
         call xerrwv(40hdresol-- neqlen too small: run aborted  ,
     1      40,1,2,0,0,0,0,0.e0,0.e0)
         return
      endif
      nns=nn*nn
      ney=neqn+2*nns
      if (ny.lt.ney) then
         istate=-3
         call xerrwv(40hdresol-- ny too small: run aborted      ,
     1      40,1,2,0,0,0,0,0.e0,0.e0)
         return
      endif
      nqs=nnq*nnq
      nps=nnp*nnp
      iwork(1)=0
      iwork(2)=0
      nflevc=1
      nflevp=1
      neqn1=neqn+1
      neqn11=neqn1+nns
      neqn12=neqn11+nqs
      neqn21=neqn12+neqn
      neqn22=neqn21+neqn
      np3=nps+3
      np23=np3+nps
      np23q=np23+nqs
      ntt=np23q+nqs
      nttp=ntt+nnp
      ntt2p=nttp+nnp
      ntt1=ntt2p+nnp
      ntt1q=ntt1+nnq
      ntt2q=ntt1q+nnq
c splitting of the various arrays completed
c ready to call the solver
 10   call drequs(neq,x,nx,y,mf,t,tout,rtol,atol,itask,istate,iopt,
     *            rwork,lrex,iwork,liw,probl,rarr,iarr)
      return
c ----------------- end of subroutine dresol -----------------------------
      end
      subroutine drequs(neq,x,nx,y,mf,t,tout,
     *                  rtol,atol,itask,istate,iopt,
     *                  rwork,lrw,iwork,liw,probl,rarr,iarr)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
      external probl,f,jac
      integer mf, itask, istate, iopt, lrw, iwork, liw, iarr
      real x, y, t, tout, rtol, atol, rwork, rarr
      dimension neq(*),x(nx,*),y(*),rwork(lrw),
     1   iwork(liw), rarr(*), iarr(*)
c subroutine drequs completes the interface to the solver drequ
c it does further array splitting and sets the input to conform
c to the LSODE solver of hindmarsh.
      integer i, j, inpj
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      real reuvct, reuvpt, tcj, tpj
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
      if (istate.ne.1.and.istate.ne.3) goto 10
      do 5 j=1,nnq
         do 5 i=1,nnp
            inpj=i+nnp*(j-1)
            y(inpj)=x(i,j)
  5   continue
c force scalar tolerances to drequ (the revised lsode)
 10   itol=1
      call drequ(f,neq,y,t,tout,itol,rtol,atol,itask,istate,iopt,
     *           rwork,lrw,iwork,liw,jac,mf,probl,rarr,iarr)
      do 14 j=1,nnq
         do 14 i=1,nnp
            inpj=i+nnp*(j-1)
 14         x(i,j)=y(inpj)
      return
c ---------- end of subroutine drequs -------------------------------
      end
      subroutine drequ (f, neq, y, t, tout, itol, rtol, atol, itask,
     1                  istate, iopt, rwork, lrw, iwork, liw, jac, mf,
     2                  probl, rarr, iarr)
c
c revision date: July 10 1990
c author of modifications: Luca Dieci, School of Mathematics, GaTech
c
c history: this is essentially "subroutine lsode" written by Alan
c          Hindmarsh at LLL, with the necessary modifications for
c          differential Riccati equations
c
      external f, jac, probl
      integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf, iarr
      real y, t, tout, rtol, atol, rwork
      dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
      dimension rarr(*), iarr(*)
c ***********************************************************
c       This program drequ constitutes a modification of the
c       lsode package made by Luca Dieci at georgia tech
c       for solving Differential Riccati Equations.
c       The modified code mantains many of the features
c       of the lsode package, but has a totally redesigned
c       nonlinear system solver and error control
c       (being based also upon the eigenvalues of the Jacobian)
c ************************************************************
c
c drequ solves the initial value problem for stiff or nonstiff
c systems of first order Riccati Differential Equations
c
c     dX/dt = A21(t) + A22(t)*X -X*A11(t) - X*A12(t)*X
c     + appropriate ICs.
c
c      The program actually rewrites (because of the
c      column orientation of fortran, but matrix arithmetic
c      is performed throughout) the riccati equation as
c
c     dy/dt = f(t,y) ,  or, in component form,
c     dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neqn)) (i = 1,...,neqn).
c
c the user, however, should never call drequ directly, but
c rather use the calling sequence as described in drequ
c the next information should be read only to use the higher capabilities
c of the dresol package, but the extra available options should not
c be necessary.
c
c-----------------------------------------------------------------------
c author and contact..
c                      (original author of the LSODE package)
c                      alan c. hindmarsh,
c                      computing and mathematics research div., l-316
c                      lawrence livermore national laboratory
c                      livermore, ca 94550.
c
c                      (author of the solver for riccati equations)
c                      luca dieci (*)
c                      school of mathematics
c                      georgia institute of technology
c                      atlanta, ga 30332
c
c (*) author to whom correspondance related to the riccati
c     equation should be addressed, and unique responsible
c     for the possible mulfunctioning of this
c     riccati solver.  please communicate promptly
c     any mistake you should find.
c
c ----------------------------------------------------------
c the descriptions of the arguments which are not examined in
c the introduction to dresol follows. we refer to the introduction
c to dresol for a full explanation of the other parameters.
c
c f      = the name of the subroutine defining the riccati ode system.
c          the system has been internally 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 through the related
c          subroutines doblks and foftex (user might want to
c          look at these routines if he desires to call drequ
c          directly).  in any case, the purpose is to provide dresol
c          with a way to compute the function f.
c	   the subroutine f has the form
c		subroutine f(t, y, ydot, probl, rarr, iarr)
c               dimension y(*), ydot(*), neq(*), rarr(*), iarr(*)
c	   where t, probl and y are input, and ydot = f(t,y) is output.
c	   for the dimensions of y needed, see later (ydot does not
c	   need to be explicitely dimensioned by the calling program as
c	   it is always made up by pieces of rwork)
c          subroutine f does not alter y(1),...,y(neq(1)).
c	   but it alters the other entries of y.
c	   f is declared external in the calling program (i.e. drequ)
c          notice that subroutine f accesses internally defined
c	   quantities in y(neq(1)+1),...
c          if the derivative dy/dt is needed, use intdy instead
c          of calling f for this purpose
c
c itol   = an indicator for the type of error control.  see
c          description below under atol.  used only for input.
c          we have forced itol=1 in the routine drequs.
c          this is reasonable for riccati equations.
c
c rtol   = a relative error tolerance parameter, either a scalar or
c          an array of length neq(1).  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(1).  input only.
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                      max-norm of ( e(i)/ewt(i) )   .le.   1,
c          where       ewt(i) = rtol(i)*abs(y(i)) + atol(i),
c          and the max-norm is max-norm(v) = max(abs(v(i), i=1,neqn)
c          aslo, 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             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          when either of these parameters is a scalar, it need not
c          be dimensioned in the user-s calling program.
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, as already explained in the
c          introduction of dresol.
c
c jac    = the name of the internally generated routine (miter=1) to
c          compute the jacobian matrix.
c          this routine sets up the block of the sylvester
c          equation to solve at each newton iteration and it has the form
c               subroutine jac (nnq, nnp, x, a011, a012, a022, icontr)
c          jac must be declared external in the calling program.
c          the routine jac is called by the subroutine prepj, and
c          we refer to the documentation in prepj for further explanation
c          of the available strategy for the newton iteration
c
c-----------------------------------------------------------------------
      external prepj, solsy
      integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
     1   mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
      integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     1   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      integer i, i1, i2, iflag, imxer, kgo, lf0,
     1   leniw, lenrw, lenwm,  mord,  mxhnl0, mxstp0, nstext
      real rowns,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      real reuvct, reuvpt, tcj, tpj
      real atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
     1   tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0,
     2   r1mach, vnorm
      dimension mord(2)
      logical ihit
c-----------------------------------------------------------------------
c the following internal common block contains
c (a) variables which are local to any subroutine but whose values must
c     be preserved between calls to the routine (own variables), and
c (b) variables which are communicated between subroutines.
c the structure of the block is as follows..  all real variables are
c listed first, followed by all integers.  within each type, the
c variables are grouped with those local to subroutine drequ first,
c then those local to subroutine stode, and finally those used
c for communication.  the block is declared in subroutines
c drequ, intdy, stode, prepj, and solsy.  groups of variables are
c replaced by dummy arrays in the common declarations in routines
c where those variables are not used.
c-----------------------------------------------------------------------
      common /ls0001/ rowns(209),
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
     3   mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
c
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
c
      data  mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
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, jump to block g and return immediately.
c-----------------------------------------------------------------------
      if (istate .lt. 1 .or. istate .gt. 3) go to 601
      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 .eq. 0) go to 603
      if (istate .eq. 2) go to 200
      go to 20
 10   init = 0
      if (tout .eq. t) go to 430
 20   ntrep = 0
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 inputs and various initializations.
c
c first check legality of the non-optional inputs neq, itol, iopt,
c mf.
c-----------------------------------------------------------------------
      if (neq(1) .le. 0) go to 604
      if (istate .eq. 1) go to 25
      if (neq(1) .gt. n) go to 605
 25   n = neq(1)
      if (itol .lt. 1 .or. itol .gt. 4) go to 606
      if (iopt .ne. 1) iopt = 0
      meth = mf/10
      miter = mf - 10*meth
      if (meth .lt. 1 .or. meth .gt. 2) go to 608
      if (miter .lt. 0 .or. miter .gt. 1) go to 608
c next process and check the optional inputs. --------------------------
      if (iopt .eq. 1) go to 40
      maxord = mord(meth)
      mxstep = mxstp0
      mxhnil = mxhnl0
      if (istate .eq. 1) h0 = 0.0e0
      hmxi = 0.0e0
      hmin = 0.0e0
      go to 60
 40   maxord = iwork(5)
      if (maxord .lt. 0) go to 611
      if (maxord .eq. 0) maxord = 100
      maxord = min0(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. 0.0e0) go to 614
 50   hmax = rwork(6)
      if (hmax .lt. 0.0e0) go to 615
      hmxi = 0.0e0
      if (hmax .gt. 0.0e0) hmxi = 1.0e0/hmax
      hmin = rwork(7)
      if (hmin .lt. 0.0e0) go to 616
c-----------------------------------------------------------------------
c set work array pointers and check lengths lrw and liw.
c pointers to segments of rwork and iwork are named by prefixing l to
c the name of the segment.  e.g., the segment yh starts at rwork(lyh).
c segments of rwork (in order) are denoted  yh, wm, ewt, savf, acor.
c-----------------------------------------------------------------------
 60   lyh = 21
      if (istate .eq. 1) nyh = n
      lwm = lyh + (maxord + 1)*nyh
      nstext=3*(nnq+nnp)+2*(nqs+nps)
      if (miter .eq. 0) lenwm = 0
      if (miter .eq. 1) lenwm = nstext + 2
      lewt = lwm + lenwm
      lsavf = lewt + n
      lacor = lsavf + n
      lenrw = lacor + n - 1
      iwork(17) = lenrw
      liwm = 1
      leniw = 20 + nnq +nnp
      if (miter .eq. 0) leniw = 20
      iwork(18) = leniw
      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. 0.0e0) go to 619
        if (atoli .lt. 0.0e0) go to 620
 70   continue
      if (istate .eq. 1) go to 100
c if istate = 3, set flag to signal parameter changes to stode. --------
      jstart = -1
      if (nq .le. maxord) go to 90
c maxord was reduced below nq.  copy yh(*,maxord+2) into savf. ---------
      do 80 i = 1,n
 80     rwork(i+lsavf-1) = rwork(i+lwm-1)
c reload wm(1) = rwork(lwm), since lwm may have changed. ---------------
 90   if (miter .gt. 0) rwork(lwm) = sqrt(uround)
      if (n .eq. nyh) go to 200
c neq was reduced.  zero part of yh to avoid undefined references. -----
      i1 = lyh + l*nyh
      i2 = lyh + (maxord + 1)*nyh - 1
      if (i1 .gt. i2) go to 200
      do 95 i = i1,i2
 95     rwork(i) = 0.0e0
      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 = r1mach(4)
      tn = t
      if (itask .ne. 4 .and. itask .ne. 5) go to 110
      tcrit = rwork(1)
      if ((tcrit - tout)*(tout - t) .lt. 0.0e0) go to 625
      if (h0 .ne. 0.0e0 .and. (t + h0 - tcrit)*h0 .gt. 0.0e0)
     1   h0 = tcrit - t
 110  jstart = 0
      if (miter .gt. 0) rwork(lwm) = sqrt(uround)
      nhnil = 0
      nst = 0
      nje = 0
      nslast = 0
      hu = 0.0e0
      nqu = 0
      ccmax = 0.3e0
      maxcor = 3
      msbp = 20
      mxncf = 10
c initial call to f.  (lf0 points to yh(*,2).) -------------------------
      lf0 = lyh + nyh
      call f (t, y, rwork(lf0), probl, rarr, iarr)
      nfe = 1
c load the initial value vector in yh. ---------------------------------
      do 115 i = 1,n
 115    rwork(i+lyh-1) = y(i)
c load and invert the ewt array.  (h is temporarily set to 1.0.) -------
      nq = 1
      h = 1.0e0
      call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
      do 120 i = 1,n
        if (rwork(i+lewt-1) .le. 0.0e0) go to 621
 120    rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
c-----------------------------------------------------------------------
c the coding below computes the step size, h0, to be attempted on the
c first step, unless the user has supplied a value for this.
c first check that tout - t differs significantly from zero.
c a scalar tolerance quantity tol is computed, as max(rtol(i))
c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted
c so as to be between 100*uround and 1.0e-3.
c then the computed value h0 is given by..
c                                      neq
c   h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2  )
c                                       1
c where   w0     = max ( abs(t), abs(tout) ),
c         f(i)   = i-th component of initial value of f,
c         ywt(i) = ewt(i)/tol  (a weight for y(i)).
c the sign of h0 is inferred from the initial values of tout and t.
c-----------------------------------------------------------------------
      if (h0 .ne. 0.0e0) go to 180
      tdist = abs(tout - t)
      w0 = amax1(abs(t),abs(tout))
      if (tdist .lt. 2.0e0*uround*w0) go to 622
      tol = rtol(1)
      if (itol .le. 2) go to 140
      do 130 i = 1,n
 130    tol = amax1(tol,rtol(i))
 140  if (tol .gt. 0.0e0) go to 160
      atoli = atol(1)
      do 150 i = 1,n
        if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
        ayi = abs(y(i))
        if (ayi .ne. 0.0e0) tol = amax1(tol,atoli/ayi)
 150    continue
 160  tol = amax1(tol,100.0e0*uround)
      tol = amin1(tol,0.001e0)
      sum = vnorm (n, rwork(lf0), rwork(lewt))
      sum = 1.0e0/(tol*w0*w0) + tol*sum**2
      h0 = 1.0e0/sqrt(sum)
      h0 = amin1(h0,tdist)
      h0 = sign(h0,tout-t)
c adjust h0 if necessary to meet hmax bound. ---------------------------
 180  rh = abs(h0)*hmxi
      if (rh .gt. 1.0e0) h0 = h0/rh
c load h with h0 and scale yh(*,2) by h0. ------------------------------
      h = h0
      do 190 i = 1,n
 190    rwork(i+lf0-1) = h0*rwork(i+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
      go to (210, 250, 220, 230, 240), itask
 210  if ((tn - tout)*h .lt. 0.0e0) go to 250
      call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
      if (iflag .ne. 0) go to 627
      t = tout
      go to 420
 220  tp = tn - hu*(1.0e0 + 100.0e0*uround)
      if ((tp - tout)*h .gt. 0.0e0) go to 623
      if ((tn - tout)*h .lt. 0.0e0) go to 250
      go to 400
 230  tcrit = rwork(1)
      if ((tn - tcrit)*h .gt. 0.0e0) go to 624
      if ((tcrit - tout)*h .lt. 0.0e0) go to 625
      if ((tn - tout)*h .lt. 0.0e0) go to 245
      call intdy (tout, 0, rwork(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. 0.0e0) go to 624
 245  hmx = abs(tn) + abs(h)
      ihit = abs(tn - tcrit) .le. 100.0e0*uround*hmx
      if (ihit) go to 400
      tnext = tn + h*(1.0e0 + 4.0e0*uround)
      if ((tnext - tcrit)*h .le. 0.0e0) go to 250
      h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
      if (istate .eq. 2) jstart = -2
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 stode.
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 ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
      do 260 i = 1,n
        if (rwork(i+lewt-1) .le. 0.0e0) go to 510
 260    rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
 270  tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt))
      if (tolsf .le. 1.0e0) go to 280
      tolsf = tolsf*2.0e0
      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
      call xerrwv(50hdresol-- warning..internal t (=r1) and h (=r2) are,
     1   50, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(
     1  60h      such that in the machine, t + h = t on the next step  ,
     1   60, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      (h = step size). solver will continue anyway,
     1   50, 101, 0, 0, 0, 0, 2, tn, h)
      if (nhnil .lt. mxhnil) go to 290
      call xerrwv(50hdresol--  above warning has been issued i1 times  ,
     1   50, 102, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      it will not be issued again for this problem,
     1   50, 102, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
 290  continue
c-----------------------------------------------------------------------
c     call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,
c    1   jac,prepj,solsy,probl,rarr,iarr)
c-----------------------------------------------------------------------
      call stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
     1   rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
     2   f, jac, prepj, solsy, probl, rarr, iarr)
      kgo = 1 - kflag
      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
      go to (310, 400, 330, 340, 350), itask
c itask = 1.  if tout has been reached, interpolate. -------------------
 310  if ((tn - tout)*h .lt. 0.0e0) go to 250
      call intdy (tout, 0, rwork(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. 0.0e0) 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. 0.0e0) go to 345
      call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
      t = tout
      go to 420
 345  hmx = abs(tn) + abs(h)
      ihit = abs(tn - tcrit) .le. 100.0e0*uround*hmx
      if (ihit) go to 400
      tnext = tn + h*(1.0e0 + 4.0e0*uround)
      if ((tnext - tcrit)*h .le. 0.0e0) go to 250
      h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
      jstart = -2
      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. 100.0e0*uround*hmx
c-----------------------------------------------------------------------
c block g.
c the following block handles all successful returns from drequ.
c if itask .ne. 1, y is loaded from yh and t is set accordingly.
c istate is set to 2, the illegal input counter is zeroed, and the
c optional outputs are loaded into the work arrays before returning.
c if istate = 1 and tout = t, there is a return with no action taken,
c except that if this has happened repeatedly, the run is terminated.
c-----------------------------------------------------------------------
 400  do 410 i = 1,n
 410    y(i) = rwork(i+lyh-1)
      t = tn
      if (itask .ne. 4 .and. itask .ne. 5) go to 420
      if (ihit) t = tcrit
 420  istate = 2
      illin = 0
      rwork(11) = hu
      rwork(12) = h
      rwork(13) = tn
      iwork(11) = nst
      iwork(12) = nfe
      iwork(13) = nje
      iwork(14) = nqu
      iwork(15) = nq
      return
c
 430  ntrep = ntrep + 1
      if (ntrep .lt. 5) return
      call xerrwv(
     1  60hdresol--  repeated calls with istate = 1 and tout = t (=r1) ,
     1   60, 301, 0, 0, 0, 0, 1, t, 0.0e0)
      go to 800
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, t is set to tn, and the illegal input
c counter illin is set to 0.  the optional outputs are loaded into
c the work arrays before returning.
c-----------------------------------------------------------------------
c the maximum number of steps was taken before reaching tout. ----------
 500  call xerrwv(50hdresol--  at current t (=r1), mxstep (=i1) steps  ,
     1   50, 201, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      taken on this call before reaching tout     ,
     1   50, 201, 0, 1, mxstep, 0, 1, tn, 0.0e0)
      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)
      call xerrwv(50hdresol-- at t (=r1), ewt(i1) has become r2 .le. 0.,
     1   50, 202, 0, 1, i, 0, 2, tn, ewti)
      istate = -6
      go to 580
c too much accuracy requested for machine precision. -------------------
 520  call xerrwv(50hdresol--  at t (=r1), too much accuracy requested ,
     1   50, 203, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      for precision of machine..  see tolsf (=r2) ,
     1   50, 203, 0, 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  call xerrwv(50hdresol-- at t(=r1) and step size h(=r2), the error,
     1   50, 204, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      test failed repeatedly or with abs(h) = hmin,
     1   50, 204, 0, 0, 0, 0, 2, tn, h)
      istate = -4
      go to 560
c kflag = -2.  convergence failed repeatedly or with abs(h) = hmin. ----
 540  call xerrwv(50hdresol--  at t (=r1) and step size h (=r2), the   ,
     1   50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(50h      corrector convergence failed repeatedly     ,
     1   50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(30h      or with abs(h) = hmin   ,
     1   30, 205, 0, 0, 0, 0, 2, tn, h)
      istate = -5
c compute imxer if relevant. -------------------------------------------
 560  big = 0.0e0
      imxer = 1
      do 570 i = 1,n
        size = abs(rwork(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, illin, and optional outputs. ------------------------
 580  do 590 i = 1,n
 590    y(i) = rwork(i+lyh-1)
      t = tn
      illin = 0
      rwork(11) = hu
      rwork(12) = h
      rwork(13) = tn
      iwork(11) = nst
      iwork(12) = nfe
      iwork(13) = nje
      iwork(14) = nqu
      iwork(15) = nq
      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.  then if there have been
c 5 consecutive such returns just before this call to the solver,
c the run is halted.
c-----------------------------------------------------------------------
 601  call xerrwv(30hdresol--  istate (=i1) illegal,
     1   30, 1, 0, 1, istate, 0, 0, 0.0e0, 0.0e0)
      go to 700
 602  call xerrwv(30hdresol--  itask (=i1) illegal ,
     1   30, 2, 0, 1, itask, 0, 0, 0.0e0, 0.0e0)
      go to 700
 603  call xerrwv(50hdresol--  istate .gt. 1 but dresol not initialized,
     1   50, 3, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      go to 700
 604  call xerrwv(30hdresol--  neq (=i1) .lt. 1    ,
     1   30, 4, 0, 1, neq(1), 0, 0, 0.0e0, 0.0e0)
      go to 700
 605  call xerrwv(50hdresol--  istate = 3 and neq increased (i1 to i2) ,
     1   50, 5, 0, 2, n, neq(1), 0, 0.0e0, 0.0e0)
      go to 700
 606  call xerrwv(30hdresol--  itol (=i1) illegal  ,
     1   30, 6, 0, 1, itol, 0, 0, 0.0e0, 0.0e0)
      go to 700
 608  call xerrwv(30hdresol--  mf (=i1) illegal    ,
     1   30, 8, 0, 1, mf, 0, 0, 0.0e0, 0.0e0)
      go to 700
 611  call xerrwv(30hdresol--  maxord (=i1) .lt. 0 ,
     1   30, 11, 0, 1, maxord, 0, 0, 0.0e0, 0.0e0)
      go to 700
 612  call xerrwv(30hdresol--  mxstep (=i1) .lt. 0 ,
     1   30, 12, 0, 1, mxstep, 0, 0, 0.0e0, 0.0e0)
      go to 700
 613  call xerrwv(30hdresol--  mxhnil (=i1) .lt. 0 ,
     1   30, 13, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
      go to 700
 614  call xerrwv(40hdresol--  tout (=r1) behind t (=r2)     ,
     1   40, 14, 0, 0, 0, 0, 2, tout, t)
      call xerrwv(50h      integration direction is given by h0 (=r1)  ,
     1   50, 14, 0, 0, 0, 0, 1, h0, 0.0e0)
      go to 700
 615  call xerrwv(30hdresol--  hmax (=r1) .lt. 0.0 ,
     1   30, 15, 0, 0, 0, 0, 1, hmax, 0.0e0)
      go to 700
 616  call xerrwv(30hdresol--  hmin (=r1) .lt. 0.0 ,
     1   30, 16, 0, 0, 0, 0, 1, hmin, 0.0e0)
      go to 700
 617  call xerrwv(
     1  60hdrequ--  rwork length needed, lenrw (=i1), exceeds lrw (=i2),
     1   60, 17, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
      go to 700
 618  call xerrwv(
     1  60hdrequ--  iwork length needed, leniw (=i1), exceeds liw (=i2),
     1   60, 18, 0, 2, leniw, liw, 0, 0.0e0, 0.0e0)
      go to 700
 619  call xerrwv(40hdrequ--  rtol(i1) is r1 .lt. 0.0        ,
     1   40, 19, 0, 1, i, 0, 1, rtoli, 0.0e0)
      go to 700
 620  call xerrwv(40hdrequ--  atol(i1) is r1 .lt. 0.0        ,
     1   40, 20, 0, 1, i, 0, 1, atoli, 0.0e0)
      go to 700
 621  ewti = rwork(lewt+i-1)
      call xerrwv(40hdresol--  ewt(i1) is r1 .le. 0.0        ,
     1   40, 21, 0, 1, i, 0, 1, ewti, 0.0e0)
      go to 700
 622  call xerrwv(
     1  60hdresol-- tout (=r1) too close to t(=r2) to start integration,
     1   60, 22, 0, 0, 0, 0, 2, tout, t)
      go to 700
 623  call xerrwv(
     1  60hdresol--  itask = i1 and tout (=r1) behind tcur - hu (= r2) ,
     1   60, 23, 0, 1, itask, 0, 2, tout, tp)
      go to 700
 624  call xerrwv(
     1  60hdresol--  itask = 4 or 5 and tcrit (=r1) behind tcur (=r2)  ,
     1   60, 24, 0, 0, 0, 0, 2, tcrit, tn)
      go to 700
 625  call xerrwv(
     1  60hdresol--  itask = 4 or 5 and tcrit (=r1) behind tout (=r2)  ,
     1   60, 25, 0, 0, 0, 0, 2, tcrit, tout)
      go to 700
 626  call xerrwv(50hdresol--  at start of problem, too much accuracy  ,
     1   50, 26, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
      call xerrwv(
     1  60h      requested for precision of machine..  see tolsf (=r1) ,
     1   60, 26, 0, 0, 0, 0, 1, tolsf, 0.0e0)
      rwork(14) = tolsf
      go to 700
 627  call xerrwv(50hdresol-- trouble from intdy. itask = i1, tout = r1,
     1   50, 27, 0, 1, itask, 0, 1, tout, 0.0e0)
c
 700  if (illin .eq. 5) go to 710
      illin = illin + 1
      istate = -3
      return
 710  call xerrwv(50hdresol--  repeated occurrences of illegal input   ,
     1   50, 302, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
c
 800  call xerrwv(50hdresol--  run aborted.. apparent infinite loop    ,
     1   50, 303, 2, 0, 0, 0, 0, 0.0e0, 0.0e0)
      return
c----------------------- end of subroutine drequ ----------------------
      end
      block data atlan1
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
c this block data is needed for initialization of the
c quantities illin and ntrep.  In the original LSODE
c these quantities were initialized via a DATA statement, but this
c is not allowed on the SUN nor on the IBM (it is fine on the CDC)
      integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
     1   mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
      integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     1   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      real rowns,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      common /ls0001/ rowns(209),
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
     3   mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      data illin /0/, ntrep /0/
c----------------------- end of block-data atlan1 ----------------------
      end
      subroutine intdy (t, k, yh, nyh, dky, iflag)
c
c author: Alan Hindmarsh, LLL, in the original LSODE package
c
      integer k, nyh, iflag
      integer iownd, iowns,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      integer i, ic, j, jb, jb2, jj, jj1, jp1
      real t, yh, dky
      real rowns,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      real c, r, s, tp
      dimension yh(nyh,*), dky(*)
      common /ls0001/ rowns(209),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     3   iownd(14), iowns(6),
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
c-----------------------------------------------------------------------
c intdy 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 = neq, 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-----------------------------------------------------------------------
      iflag = 0
      if (k .lt. 0 .or. k .gt. nq) go to 80
      tp = tn - hu -  100.0e0*uround*(tn + hu)
      if ((t-tp)*(t-tn) .gt. 0.0e0) 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 = float(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 = float(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)
      do 60 i = 1,n
 60     dky(i) = r*dky(i)
      return
c
 80   call xerrwv(30hintdy--  k (=i1) illegal      ,
     1   30, 51, 0, 1, k, 0, 0, 0.0e0, 0.0e0)
      iflag = -1
      return
 90   call xerrwv(30hintdy--  t (=r1) illegal      ,
     1   30, 52, 0, 0, 0, 0, 1, t, 0.0e0)
      call xerrwv(
     1  60h      t not in interval tcur - hu (= r1) to tcur (=r2)      ,
     1   60, 52, 0, 0, 0, 0, 2, tp, tn)
      iflag = -2
      return
c----------------------- end of subroutine intdy -----------------------
      end
      subroutine stode (neq, y, yh, nyh, yh1, ewt, savf, acor,
     1   wm, iwm, f, jac, pjac, slvs, probl, rarr, iarr)
c
c revision date: July 10 1990
c modifications' author: Luca Dieci, School of Mathematics, GaTech
c
c history: based on the routine by the same name, which is in the
c          lsode package of A.Hindmarsh.  This version has many
c          changes
c
      external f, jac, pjac, slvs, probl
      integer neq, nyh, iwm, iarr
      integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      integer i, i1, iredo, iret, j, jb, m, ncf, newq
      integer ksup
      real y, yh, yh1, ewt, savf, acor, wm, rarr
      real conit, crate, el, elco, hold, rmax, tesco,
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      real reuvct, reuvpt, tcj, tpj
      real dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
     1   r, rh, rhdn, rhsm, rhup, told, vnorm
      real tpjph, signh
      dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
     1   acor(*), wm(*), rarr(*), iwm(*), iarr(*)
      common /ls0001/ conit, crate, el(13), elco(13,12),
     1   hold, rmax, tesco(3,12),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14),
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
c-----------------------------------------------------------------------
c stode performs one step of the integration of an initial value
c problem for a system of ordinary differential equations.
c note.. stode is independent of the value of the iteration method
c indicator miter, when this is .ne. 0, and hence is independent
c of the type of chord method used.
c communication with stode is done with the following variables..
c
c neq    = integer array containing problem size in neq(1), and
c          passed as the neq argument in all calls to f and jac.
c y      = an array of length .ge. n used as the y argument in
c          all calls to f and jac.
c yh     = an nyh 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 nyh    = a constant integer .ge. n, the first dimension of yh.
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 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 = real and integer work arrays associated with matrix
c          operations in chord iteration (miter .ne. 0).
c pjac   = name of routine to evaluate and preprocess jacobian matrix
c          and p = i - h*el0*jac, if a chord method is being used.
c slvs   = name of routine to solve linear system in chord iteration.
c    N.B.: the routines pjac and slvs actually consist in simple
c          interfaces calling various other routines.  See their
c          own documentation
c ccmax  = maximum relative change in h*el0 before pjac is called.
c h      = the step size to be attempted on the next step.
c          h is altered by the error control algorithm during the
c          problem.  h can be either positive or negative, but its
c          sign must remain constant throughout the problem.
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          hmin and hmxi may be changed at any time, but will not
c          take effect until the next change of h is considered.
c tn     = the independent variable. tn is updated on each step taken.
c jstart = an integer used for input only, with the following
c          values and meanings..
c               0  perform the first step.
c           .gt.0  take a new step continuing from the last.
c              -1  take the next step with a new value of h, maxord,
c                    n, meth, miter, and/or matrix parameters.
c              -2  take the next step with a new value of h,
c                    but with other inputs unchanged.
c          on return, jstart is set to 1 to facilitate continuation.
c kflag  = a completion code 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  fatal error in pjac or slvs.
c          a return with kflag = -1 or -2 means either
c          abs(h) = hmin or 10 consecutive failures occurred.
c          on a return with kflag negative, the values of tn and
c          the yh array are as of the beginning of the last
c          step, and h is the last step size attempted.
c maxord = the maximum order of integration method to be allowed.
c maxcor = the maximum number of corrector iterations allowed.
c msbp   = maximum number of steps between pjac calls (miter .gt. 0).
c mxncf  = maximum number of convergence failures allowed.
c meth/miter = the method flags.  see description in driver.
c n      = the number of first-order differential equations.
c
c --- a number of additions have been made in this routine with respect
c  the original one of LSODE.  In particular, we handle here the option
c  of how to solve the nolinear system and the superstability option/check
c  based upon eigenvalue information.  Communication is provided by the
c  following variables in the block /intdre/
c
c reuvct = most unstable real part of the eigenvalues at current time t
c reuvpt = most unstable real part of the eigenvalues at previous time t
c tcj = time at which the most current evaluation/factorization of the
c       Jacobian matrix has been made
c tpj = next to the last time at which an evaluation/factorization of the
c       Jacobian martix has been made
c nsuper = integer flag signalling whether the code should (nsuper = 1)
c          perform a superstability check or not (nsuper = 0).  For the
c          criterion adopted we refer to the documentation of DRESOL
c nflevc = integer flaf monitoring whether at the current step we had
c          (nflevc = 1) or not (nflevc = 0) unstable eigenvalues with
c          respect to the direction of integration and the machine precision
c nflevp = integer flag as nflevc, but recording unstable/stable eigenvalues
c          at the previous step
c nonlin = the way we should solve the nonlinear system of equations in
c          case miter = 1 (i.e. mf=11 or 21).  We set nonlin = 1 if we
c          want a real quasi-Newton iteration (that is, a factorizatio
c          at each step and then freeze the jacobian for that one step)
c          and we set nonlin = 0 if we want a chord iteration with a
c          factored jacobian being kept across many steps (depending on msbp)
c
c all of the other parameters are as explained in dresol
c-----------------------------------------------------------------------
      kflag = 0
      told = tn
      ncf = 0
      ierpj = 0
      iersl = 0
      jcur = 0
      icf = 0
      delp = 0.0e0
      if (jstart .eq. 0) tcj = tn
      tpj = tcj
      nflevp = nflevc
      reuvpt = reuvct
      ksup = 0
      if (jstart .gt. 0) go to 200
      if (jstart .eq. -1) go to 100
      if (jstart .eq. -2) go to 160
c-----------------------------------------------------------------------
c on the first call, the order is set to 1, and other variables are
c initialized.  rmax is the maximum ratio by which h can be increased
c in a single step.  it is initially 1.e4 to compensate for the small
c initial h, but then is normally equal to 10.  if a failure
c occurs (in corrector convergence or error test), rmax is set at 2
c for the next increase.
c-----------------------------------------------------------------------
      lmax = maxord + 1
      nq = 1
      l = 2
      ialth = 2
      rmax = 10000.0e0
      rc = 0.0e0
      el0 = 1.0e0
      crate = 0.7e0
      hold = h
      meo = meth
      nslp = 0
      ipup = miter
      iret = 3
      go to 140
c-----------------------------------------------------------------------
c the following block handles preliminaries needed when jstart = -1.
c ipup is set to miter to force a matrix update.
c if an order increase is about to be considered (ialth = 1),
c ialth is reset to 2 to postpone consideration one more step.
c if the caller has changed meth, cfode is called to reset
c the coefficients of the method.
c if the caller has changed maxord to a value less than the current
c order nq, nq is reduced to maxord, and a new h chosen accordingly.
c if h is to be changed, yh must be rescaled.
c if h or meth is being changed, ialth is reset to l = nq + 1
c to prevent further changes in h for that many steps.
c-----------------------------------------------------------------------
 100  ipup = miter
      lmax = maxord + 1
      if (ialth .eq. 1) ialth = 2
      if (meth .eq. meo) go to 110
      call cfode (meth, elco, tesco)
      meo = meth
      if (nq .gt. maxord) go to 120
      ialth = l
      iret = 1
      go to 150
 110  if (nq .le. maxord) go to 160
 120  nq = maxord
      l = lmax
      do 125 i = 1,l
 125    el(i) = elco(i,nq)
      nqnyh = nq*nyh
      rc = rc*el(1)/el0
      el0 = el(1)
      conit = 0.5e0/float(nq+2)
      ddn = vnorm (n, savf, ewt)/tesco(1,l)
      exdn = 1.0e0/float(l)
      rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
      rh = amin1(rhdn,1.0e0)
      iredo = 3
      if (h .eq. hold) go to 170
      rh = amin1(rh,abs(h/hold))
      h = hold
      go to 175
c-----------------------------------------------------------------------
c cfode is called to get all the integration coefficients for the
c current meth.  then the el vector and related constants are reset
c whenever the order nq is changed, or at the start of the problem.
c-----------------------------------------------------------------------
 140  call cfode (meth, elco, tesco)
 150  do 155 i = 1,l
 155    el(i) = elco(i,nq)
      nqnyh = nq*nyh
      rc = rc*el(1)/el0
      el0 = el(1)
      conit = 0.5e0/float(nq+2)
      go to (160, 170, 200), iret
c-----------------------------------------------------------------------
c if h is being changed, the h ratio rh is checked against
c rmax, hmin, and hmxi, and the yh array rescaled.  ialth is set to
c l = nq + 1 to prevent a change of h for that many steps, unless
c forced by a convergence or error test failure.
c-----------------------------------------------------------------------
 160  if (h .eq. hold) go to 200
      rh = h/hold
      h = hold
      iredo = 3
      go to 175
 170  rh = amax1(rh,hmin/abs(h))
 175  rh = amin1(rh,rmax)
      rh = rh/amax1(1.0e0,abs(h)*hmxi*rh)
      r = 1.0e0
      do 180 j = 2,l
        r = r*rh
        do 180 i = 1,n
 180      yh(i,j) = yh(i,j)*r
      h = h*rh
      rc = rc*rh
      ialth = l
      if (iredo .eq. 0) go to 690
c-----------------------------------------------------------------------
c this section computes the predicted values by effectively
c multiplying the yh array by the pascal triangle matrix.
c rc is the ratio of new to old values of the coefficient  h*el(1).
c when rc differs from 1 by more than ccmax, ipup is set to miter
c to force pjac to be called, if a jacobian is involved.
c in any case, pjac is called at least every msbp steps.
c if nonlin = 1 pjac is called every step.
c-----------------------------------------------------------------------
 200  if (abs(rc-1.0e0) .gt. ccmax) ipup = miter
      if (nst .ge. nslp+msbp) ipup = miter
c next line forces a genuine quasi-newton iteration, with
c the jacobian updated at each new step, if nonlin=1
      if (nonlin .eq. 1) ipup = miter
      tn = tn + h
      i1 = nqnyh + 1
      do 215 jb = 1,nq
        i1 = i1 - nyh
        do 210 i = i1,nqnyh
 210      yh1(i) = yh1(i) + yh1(i+nyh)
 215    continue
c-----------------------------------------------------------------------
c up to maxcor corrector iterations are taken.  a convergence test is
c made on the max norm of each correction, weighted by the error
c weight vector ewt.  the sum of the corrections is accumulated in the
c vector acor(i).  the yh array is not altered in the corrector loop.
c-----------------------------------------------------------------------
 220  m = 0
      do 230 i = 1,n
 230    y(i) = yh(i,1)
      call f (tn, y, savf, probl, rarr, iarr)
      nfe = nfe + 1
      if (ipup .le. 0) go to 250
c-----------------------------------------------------------------------
c if indicated, the matrix p = i - h*el(1)*j is reevaluated and
c preprocessed before starting the corrector iteration.  ipup is set
c to 0 as an indicator that this has been done.  see pjac for how this is done.
c-----------------------------------------------------------------------
      call pjac (y,wm,iwm,probl,jac,rarr,iarr)
      ipup = 0
      rc = 1.0e0
      nslp = nst
      crate = 0.7e0
      if (ierpj .ne. 0) go to 430
 250  do 260 i = 1,n
 260    acor(i) = 0.0e0
 270  if (miter .ne. 0) go to 350
c-----------------------------------------------------------------------
c in the case of functional iteration, update y directly from
c the result of the last function evaluation.
c-----------------------------------------------------------------------
      do 290 i = 1,n
        savf(i) = h*savf(i) - yh(i,2)
 290    y(i) = savf(i) - acor(i)
      del = vnorm (n, y, ewt)
      do 300 i = 1,n
        y(i) = yh(i,1) + el(1)*savf(i)
 300    acor(i) = savf(i)
      go to 400
c-----------------------------------------------------------------------
c in the case of the chord or quasi-Newton method, compute the corrector
c error, and solve the linear system with that as right-hand side and
c p as coefficient matrix.
c-----------------------------------------------------------------------
 350  do 360 i = 1,n
 360    y(i) = h*savf(i) - (yh(i,2) + acor(i))
      call slvs (wm, y)
      if (iersl .lt. 0) go to 430
      if (iersl .gt. 0) go to 410
      del = vnorm (n, y, ewt)
      do 380 i = 1,n
        acor(i) = acor(i) + y(i)
 380    y(i) = yh(i,1) + el(1)*acor(i)
c-----------------------------------------------------------------------
c test for convergence.  if m.gt.0, an estimate of the convergence
c rate constant is stored in crate, and this is used in the test.
c-----------------------------------------------------------------------
 400  continue
      if (m .ne. 0) crate = amax1(0.2e0*crate,del/delp)
      dcon = del*amin1(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
      if (dcon .le. 1.0e0) go to 450
      m = m + 1
      if (m .eq. maxcor) go to 410
      if (m .ge. 2 .and. del .gt. 2.0e0*delp) go to 410
      delp = del
      call f (tn, y, savf, probl, rarr, iarr)
      nfe = nfe + 1
      go to 270
c-----------------------------------------------------------------------
c the corrector iteration failed to converge.
c if miter .ne. 0 and the jacobian is out of date, pjac is called for
c the next try.  otherwise the yh array is retracted to its values
c before prediction, and h is reduced, if possible.  if h cannot be
c reduced or mxncf failures have occurred, exit with kflag = -2.
c-----------------------------------------------------------------------
 410  if (miter .eq. 0 .or. jcur .eq. 1) go to 430
      icf = 1
      ipup = miter
      go to 220
 430  icf = 2
      ncf = ncf + 1
      rmax = 2.0e0
      tn = told
      i1 = nqnyh + 1
      do 445 jb = 1,nq
        i1 = i1 - nyh
        do 440 i = i1,nqnyh
 440      yh1(i) = yh1(i) - yh1(i+nyh)
 445    continue
      if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680
      if (abs(h) .le. hmin*1.00001e0) go to 670
      if (ncf .eq. mxncf) go to 670
      rh = 0.25e0
      ipup = miter
      iredo = 1
      go to 170
c-----------------------------------------------------------------------
c the corrector has converged.  jcur is set to 0
c to signal that the jacobian involved may need updating later.
c the local error test is made and control passes to statement 500
c if it fails.
c-----------------------------------------------------------------------
 450  jcur = 0
      if (m .eq. 0) dsm = del/tesco(2,nq)
      if (m .gt. 0) dsm = vnorm (n, acor, ewt)/tesco(2,nq)
      if (dsm .gt. 1.0e0) go to 500
c-----------------------------------------------------------------------
c after a successful step, update the yh array.
c consider changing h if ialth = 1.  otherwise decrease ialth by 1.
c if ialth is then 1 and nq .lt. maxord, then acor is saved for
c use in a possible order increase on the next step.
c if a change in h is considered, an increase or decrease in order
c by one is considered also.  a change in h is made only if it is by a
c factor of at least 1.1.  if not, ialth is set to 3 to prevent
c testing for that many steps.
c-----------------------------------------------------------------------
c	       ----- superstability/eigenvalues check -----
c the last decomposition of the jacobian did not occur at current t,
c or we do not have the eigenvalues altogether
      if (tn .ne. tcj .or. miter .eq. 0) goto 466
c ready to check the sign of eigenvalues of jacobian
c do this regardless of superstability check, to signal
c possible instabilities in integration to the user
c in particular, iwm(1) = 1 means purely imaginary e.vs.  it is not
c error but user might want to use adams formulas if he is using
c the bdf.  iwm(2) = 1 means unstable real parts of eigenvalues
c with respect to direction of integration: sign(h) * real-part(e.vs) > 0
c initially, have iwm(1) = iwm(2) = 0.
      if (iwm(1) .eq. 1) goto 457
      iwm(1) = 20
      if (icontr .eq. 1) goto 454
      do 452 i = 1,nnp
            if ((iwm(20+i) .ne. 0) .and. (wm(nttp+i-1) .eq. 0))
     1           iwm(1) = 10
 452  continue
      goto 457
 454  do 455 i = 1,nnp
         do 455 j = 1,nnq
            if ((iwm(20+i) + iwm(20+nnp+j) .ne. 0) .and.
     1   (wm(nttp+i-1) - wm(ntt1q+j-1) .eq. 0)) iwm(1) = 10
 455  continue
c above, we were looking at itypea(i), itypeb(j) of the routine hqr3 (see pjac)
 457  tpjph = tpj + h
      if ((iwm(2) .eq. 1) .and. (nsuper .eq. 0 .or.
     1   tcj .ne. tpjph)) goto 466
      if (iwm(2) .ne. 1) iwm(2) = 20
      signh = -1.e0
      if (h .ge. 0.e0) signh = 1.e0
      nflevc = 0
      if (signh .lt. 0.e0) goto 460
c we check whether the most unstable eigenvalue is unstable with
c respect to the direction of integration
c in case of icontr=0, the eigenvalues are ordered in incresing
c order of real parts, and they have the reversed sign.
      if (icontr .eq. 0) reuvct = - wm(nttp)
      if (icontr .eq. 1) reuvct = wm(nttp+nnp-1)-wm(ntt1q)
      goto 464
 460  if (icontr .eq. 0) reuvct = - wm(nttp+nnp-1)
      if (icontr .eq. 1) reuvct = wm(nttp)-wm(ntt1q+nnq-1)
 464  reuvct = signh * reuvct
      if (reuvct .gt. 1.0e-3) nflevc = 1
      if (nflevc .eq. 1 .and. iwm(2) .ne. 1) iwm(2) = 10
      if (nsuper .eq. 0 .or. tcj .ne. tpjph) goto 466
c check for possible superstable behavior if using the BDF or Adams
c with the Newton-like iteration (no e.vs are available, when
c using the functional fixed point iteration).
c the logic behind this check relies on (i) the assumed existence of
c a stable trajectory, or (ii) that, if we are along an unstable trajectory,
c then we have to proceed cautiously; in particular, we are not allowed
c to increase the stepsize if the degree of instability (as measured
c by the most unstable real part of the eigenvalues) is increased
c by 100% or more. If neither of these two things is true, the code might
c end up working hard for no reason at all.  We have judged as unstable
c an eigenvalue if it is positive w.r.t. h, within heuristic accuracy
c bound (1.0e-3, above).
      if (nflevc .eq. 1 .and. (abs(h) .ge. abs(hu) .or.
     1   ksup .gt. 0) .and. ((nflevp .eq. 0 .and. reuvct .gt. 1.e-2)
     2   .or. reuvct .gt. (2.0e0*reuvpt))) then
         ksup = ksup+1
         goto 465
      endif
c in the process of avoiding superstability, we have reduced the stepsize
c within machine precision without recovering from the possible occurrence
c of the phenomenon. To avoid further testing, we set nsuper to 0.  The
c user should monitor this, if wants to change it again to snuper=1 later
      if (tpjph .eq. tpj) then
         call xerrwv(47hdresol-- in subroutine stode reset nsuper to 0 ,
     1      47,11,1,0,0,0,0,0.e0,0.e0)
         call xerrwv(47h    as after i1 reductions h (=r1) is too small,
     1      47,12,1,1,ksup,0,1,h,0.e0)
         nsuper = 0
         neq(6) = 0
      endif
      ksup = 0
c if cannot recover from superstable behavior, after 6 steps reductions
c signal it to the user by letting nsuper = 0
 465  if (ksup .gt. 6) then
         call xerrwv(47hdresol-- in subroutine stode reset nsuper to 0 ,
     1      47,13,1,0,0,0,0,0.e0,0.e0)
         call xerrwv(47h    after 6 tries still superstability possible,
     1      47,14,1,0,0,0,0,0.e0,0.e0)
         nsuper = 0
         neq(6) = 0
         ksup = 0
      endif
      if (ksup .eq. 0) goto 466
c above, either the superstability check failed, or we made it more
c than 6 times without recovering.  in these cases, accept the step
c because the error tols were fine.  otherwise, force a step-size
c reduction.  after 3 failures (because of ksup or err-tols) we force
c a reduction by at least a factor of 10 at order 1.
c if we detect superstability, the step-size reduction is enforced  by
c setting dsm=10.  this should give -at first- reduction at same order
      dsm = 10
      goto 500
c ----- end of superstability/eigenvalues' check -----
 466  kflag = 0
      iredo = 0
      nst = nst + 1
      hu = h
      nqu = nq
      do 470 j = 1,l
        do 470 i = 1,n
 470      yh(i,j) = yh(i,j) + el(j)*acor(i)
      ialth = ialth - 1
      if (ialth .eq. 0) go to 520
      if (ialth .gt. 1) go to 700
      if (l .eq. lmax) go to 700
      do 490 i = 1,n
 490    yh(i,lmax) = acor(i)
      go to 700
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 this or
c one lower order.  after 2 or more failures, h is forced to decrease
c by a factor of 0.2 or less.
c-----------------------------------------------------------------------
 500  kflag = kflag - 1
      tn = told
      i1 = nqnyh + 1
      do 515 jb = 1,nq
        i1 = i1 - nyh
        do 510 i = i1,nqnyh
 510      yh1(i) = yh1(i) - yh1(i+nyh)
 515    continue
      rmax = 2.0e0
      if (abs(h) .le. hmin*1.00001e0) go to 660
      if (kflag .le. -3) go to 640
      iredo = 2
      rhup = 0.0e0
      go to 540
c-----------------------------------------------------------------------
c regardless of the success or failure of the step, factors
c rhdn, rhsm, and rhup are computed, by which h could be multiplied
c at order nq - 1, order nq, or order nq + 1, respectively.
c in the case of failure, rhup = 0.0 to avoid an order increase.
c the largest of these is determined and the new order chosen
c accordingly.  if the order is to be increased, we compute one
c additional scaled derivative.
c-----------------------------------------------------------------------
 520  rhup = 0.0e0
      if (l .eq. lmax) go to 540
      do 530 i = 1,n
 530    savf(i) = acor(i) - yh(i,lmax)
      dup = vnorm (n, savf, ewt)/tesco(3,nq)
      exup = 1.0e0/float(l+1)
      rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
 540  exsm = 1.0e0/float(l)
      rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
      rhdn = 0.0e0
      if (nq .eq. 1) go to 560
      ddn = vnorm (n, yh(1,l), ewt)/tesco(1,nq)
      exdn = 1.0e0/float(nq)
      rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
 560  if (rhsm .ge. rhup) go to 570
      if (rhup .gt. rhdn) go to 590
      go to 580
 570  if (rhsm .lt. rhdn) go to 580
      newq = nq
      rh = rhsm
      go to 620
 580  newq = nq - 1
      rh = rhdn
      if (kflag .lt. 0 .and. rh .gt. 1.0e0) rh = 1.0e0
      go to 620
 590  newq = l
      rh = rhup
      if (rh .lt. 1.1e0) go to 610
      r = el(l)/float(l)
      do 600 i = 1,n
 600    yh(i,newq+1) = acor(i)*r
      go to 630
 610  ialth = 3
      go to 700
 620  if ((kflag .eq. 0) .and. (rh .lt. 1.1e0)) go to 610
      if (kflag .le. -2) rh = amin1(rh,0.2e0)
c-----------------------------------------------------------------------
c if there is a change of order, reset nq, l, and the coefficients.
c in any case h is reset according to rh and the yh array is rescaled.
c then exit from 690 if the step was ok, or redo the step otherwise.
c-----------------------------------------------------------------------
      if (newq .eq. nq) go to 170
 630  nq = newq
      l = nq + 1
      iret = 2
      go to 150
c-----------------------------------------------------------------------
c control reaches this section if 3 or more failures have occured.
c if 10 failures have occurred, exit with kflag = -1.
c it is assumed that the derivatives that have accumulated in the
c yh array have errors of the wrong order.  hence the first
c derivative is recomputed, and the order is set to 1.  then
c h is reduced by a factor of 10, and the step is retried,
c until it succeeds or h reaches hmin.
c-----------------------------------------------------------------------
 640  if (kflag .eq. -10) go to 660
      rh = 0.1e0
      rh = amax1(hmin/abs(h),rh)
      h = h*rh
      do 645 i = 1,n
 645    y(i) = yh(i,1)
      call f (tn, y, savf, probl, rarr, iarr)
      nfe = nfe + 1
      do 650 i = 1,n
 650    yh(i,2) = h*savf(i)
      ipup = miter
      ialth = 5
      if (nq .eq. 1) go to 200
      nq = 1
      l = 2
      iret = 3
      go to 150
c-----------------------------------------------------------------------
c all returns are made through this section.  h is saved in hold
c to allow the caller to change h on the next step.
c-----------------------------------------------------------------------
 660  kflag = -1
      go to 720
 670  kflag = -2
      go to 720
 680  kflag = -3
      go to 720
 690  rmax = 10.0e0
 700  r = 1.0e0/tesco(2,nqu)
      do 710 i = 1,n
 710    acor(i) = acor(i)*r
 720  hold = h
      jstart = 1
      if (iwm(1) .eq. 1 .and. iwm(2) .eq. 1) return
      if (iwm(1) .eq. 10) iwm(1) = 1
      if (iwm(1) .eq. 20) iwm(1) = 0
      if (iwm(2) .eq. 10) iwm(2) = 1
      if (iwm(2) .eq. 20) iwm(2) = 0
      return
c----------------------- end of subroutine stode -----------------------
      end
      subroutine cfode (meth, elco, tesco)
c
c author: Alan Hindmarsh, LLL
c
      integer meth
      integer i, ib, nq, nqm1, nqp1
      real elco, tesco
      real agamq, fnq, fnqm1, pc, pint, ragq,
     1   rqfac, rq1fac, tsign, xpin
      dimension elco(13,12), tesco(3,12)
c-----------------------------------------------------------------------
c cfode is called by the integrator routine to set coefficients
c needed there.  the coefficients for the current method, as
c given by the value of meth, are set for all orders and saved.
c the maximum order assumed here is 12 if meth = 1 and 5 if meth = 2.
c (a smaller value of the maximum order is also allowed.)
c cfode is called once at the beginning of the problem,
c and is not called again unless and until meth is changed.
c
c the elco array contains the basic method coefficients.
c the coefficients el(i), 1 .le. i .le. nq+1, for the method of
c order nq are stored in elco(i,nq).  they are given by a genetrating
c polynomial, i.e.,
c     l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
c for the implicit adams methods, l(x) is given by
c     dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1),    l(-1) = 0.
c for the bdf methods, l(x) is given by
c     l(x) = (x+1)*(x+2)* ... *(x+nq)/k,
c where         k = factorial(nq)*(1 + 1/2 + ... + 1/nq).
c
c the tesco array contains test constants used for the
c local error test and the selection of step size and/or order.
c at order nq, tesco(k,nq) is used for the selection of step
c size at order nq - 1 if k = 1, at order nq if k = 2, and at order
c nq + 1 if k = 3.
c-----------------------------------------------------------------------
      dimension pc(12)
c
      go to (100, 200), meth
c
 100  elco(1,1) = 1.0e0
      elco(2,1) = 1.0e0
      tesco(1,1) = 0.0e0
      tesco(2,1) = 2.0e0
      tesco(1,2) = 1.0e0
      tesco(3,12) = 0.0e0
      pc(1) = 1.0e0
      rqfac = 1.0e0
      do 140 nq = 2,12
c-----------------------------------------------------------------------
c the pc array will contain the coefficients of the polynomial
c     p(x) = (x+1)*(x+2)*...*(x+nq-1).
c initially, p(x) = 1.
c-----------------------------------------------------------------------
        rq1fac = rqfac
        rqfac = rqfac/float(nq)
        nqm1 = nq - 1
        fnqm1 = float(nqm1)
        nqp1 = nq + 1
c form coefficients of p(x)*(x+nq-1). ----------------------------------
        pc(nq) = 0.0e0
        do 110 ib = 1,nqm1
          i = nqp1 - ib
 110      pc(i) = pc(i-1) + fnqm1*pc(i)
        pc(1) = fnqm1*pc(1)
c compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
        pint = pc(1)
        xpin = pc(1)/2.0e0
        tsign = 1.0e0
        do 120 i = 2,nq
          tsign = -tsign
          pint = pint + tsign*pc(i)/float(i)
 120      xpin = xpin + tsign*pc(i)/float(i+1)
c store coefficients in elco and tesco. --------------------------------
        elco(1,nq) = pint*rq1fac
        elco(2,nq) = 1.0e0
        do 130 i = 2,nq
 130      elco(i+1,nq) = rq1fac*pc(i)/float(i)
        agamq = rqfac*xpin
        ragq = 1.0e0/agamq
        tesco(2,nq) = ragq
        if (nq .lt. 12) tesco(1,nqp1) = ragq*rqfac/float(nqp1)
        tesco(3,nqm1) = ragq
 140    continue
      return
c
 200  pc(1) = 1.0e0
      rq1fac = 1.0e0
      do 230 nq = 1,5
c-----------------------------------------------------------------------
c the pc array will contain the coefficients of the polynomial
c     p(x) = (x+1)*(x+2)*...*(x+nq).
c initially, p(x) = 1.
c-----------------------------------------------------------------------
        fnq = float(nq)
        nqp1 = nq + 1
c form coefficients of p(x)*(x+nq). ------------------------------------
        pc(nqp1) = 0.0e0
        do 210 ib = 1,nq
          i = nq + 2 - ib
 210      pc(i) = pc(i-1) + fnq*pc(i)
        pc(1) = fnq*pc(1)
c store coefficients in elco and tesco. --------------------------------
        do 220 i = 1,nqp1
 220      elco(i,nq) = pc(i)/pc(2)
        elco(2,nq) = 1.0e0
        tesco(1,nq) = rq1fac
        tesco(2,nq) = float(nqp1)/elco(1,nq)
        tesco(3,nq) = float(nq+2)/elco(1,nq)
        rq1fac = rq1fac/fnq
 230    continue
      return
c----------------------- end of subroutine cfode -----------------------
      end
c
      subroutine f(t,y,ydot,probl,rarr,iarr)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
      external probl
      integer iarr
      real t, y, ydot, rarr
      dimension y(*),ydot(*),rarr(*),iarr(*)
c
c we give the function here
c
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      real reuvct, reuvpt, tcj, tpj
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
      call probl(t,y(neqn1),rarr,iarr)
      call doblks(nn,nnq,nnp,y(neqn1),y(neqn11),y(neqn12),
     *            y(neqn21),y(neqn22))
      call foftex(nn,nnq,nnp,y(neqn11),y(neqn12),y(neqn21),y(neqn22),
     *            y(1),ydot,y(neqn1),icontr)
      return
      end
c
      subroutine doblks(nn,nnq,nnp,a,a11,a12,a21,a22)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
      integer nn, nnq, nnp
      real a, a11, a12, a21, a22
      dimension a(nn,nn),a11(nnq,nnq),a12(nnq,nnp),
     *          a21(nnp,nnq),a22(nnp,nnp)
c
c form blocks of matrix here
c
      integer i, j
      do 5 j=1,nnq
         do 5 i=1,nnq
 5          a11(i,j)=a(i,j)
      do 10 j=1,nnp
         do 10 i=1,nnq
 10         a12(i,j)=a(i,nnq+j)
      do 15 j=1,nnq
         do 15 i=1,nnp
 15         a21(i,j)=a(nnq+i,j)
      do 20 j=1,nnp
         do 20 i=1,nnp
 20         a22(i,j)=a(nnq+i,nnq+j)
      return
      end
c
      subroutine foftex(n,nnq,nnp,a11,a12,a21,a22,
     1       xsave,effe,como,icontr)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
      integer n,nnq,nnp,icontr
      real a11, a12, a21, a22, xsave, effe, como
      dimension a11(nnq,nnq),a12(nnq,nnp),a21(nnp,nnq),a22(nnp,nnp),
     *   xsave(nnp,nnq),como(n,n),effe(nnp,nnq)
c here we form F(t,X)=A21+A22*X-X*A11-X*A12*X, where all quantities have
c been evaluated at t beforehand.
c we attempt to economize on forming F(t,X) according to whether
c NP.GT.NQ or not, and if ICONTR=0,1
      integer i, j, k
      real prec
c
      if (nnp.lt.nnq) goto 50
      do 10 j=1,nnq
         do 10 i=1,nnq
            prec=0.e0
            do 5 k=1,nnp
               if (icontr.eq.1) prec=prec+a12(i,k)*xsave(k,j)
               if (icontr.eq.0) prec=prec+a12(i,k)*xsave(k,j)*.5
 5          continue
            prec=prec+a11(i,j)
            como(i,j)=prec
 10   continue
      do 20 j=1,nnq
         do 20 i=1,nnp
            prec=0.e0
            do 15 k=1,nnq
               prec=prec-xsave(i,k)*como(k,j)
 15         continue
            effe(i,j)=prec
 20   continue
      if (icontr.eq.0) goto 35
      do 30 j=1,nnq
         do 30 i=1,nnp
            prec=0.e0
            do 25 k=1,nnp
               prec=prec+a22(i,k)*xsave(k,j)
 25         continue
            effe(i,j)=effe(i,j)+prec+a21(i,j)
 30   continue
      return
 35   continue
      do 40 j=1,nnq
         do 40 i=1,nnp
            como(i,j)=effe(j,i)
 40   continue
      do 45 j=1,nnq
         do 45 i=1,nnp
            effe(i,j)=effe(i,j)+como(i,j)+a21(i,j)
 45   continue
      return
 50   continue
      do 60 j=1,nnp
         do 60 i=1,nnp
            prec=0.e0
            do 55 k=1,nnq
               prec=prec+xsave(i,k)*a12(k,j)
 55         continue
            prec=a22(i,j)-prec
            como(i,j)=prec
 60   continue
      do 70 j=1,nnq
         do 70 i=1,nnp
            prec=0.e0
            prec=a21(i,j)
            do 65 k=1,nnp
               prec=prec+como(i,k)*xsave(k,j)
 65         continue
            effe(i,j)=prec
 70   continue
      do 80 j=1,nnq
         do 80 i=1,nnp
            prec=0.e0
            do 75 k=1,nnq
               prec=prec+xsave(i,k)*a11(k,j)
 75         continue
            effe(i,j)=effe(i,j)-prec
 80   continue
      return
      end
c
      subroutine prepj (y, wm, iwm, probl, jac, rarr, iarr)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
c history: a routine by the same name is in the original lsode of
c          Alan Hindmarsh (LLL), but this is basically all new
c
      external probl, jac
      integer iwm, iarr
      integer iownd, iowns,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      integer i, ier,  j,  lenp, nnqp1, npp1
      real y, wm, rarr
      real rowns,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      real reuvct, reuvpt, tcj, tpj
      real con, hl0
      dimension y(*), wm(*), rarr(*), iwm(*), iarr(*)
      common /ls0001/ rowns(209),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     3   iownd(14), iowns(6),
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
c-----------------------------------------------------------------------
c prepj is called by stode to compute and process the matrix
c p = i - h*el(1)*j , where j refers to the jacobian.
c the jacobian j does not need to be supplied, as it
c is (exactly) computed from the Riccati equation itself,
c via the routine jac
c
c the computation of the jacobian is separated in two distinct
c phases: (i) the computation of the blocks of the Sylvester equation
c ~         ~
c A22*X - X*A11 = C , to solve, and their decomposition , and
c (ii) the update of these blocks by premultiplication by h*el(1).
c (i) this phase is done by first forming the original matrix blocks
c (with calls to probl and doblks) and then forming the blocks of
c the transformed matrix as A22-X*A12 and A11+A12*X respectively,
c where X is the approximate solution as computed by the predictor
c formulas. The decomposition is done by our version of the
c Bartels-Stewart algorithm (with calls to axpxbd or atxxad),
c so that the transformed blocks are eigendecomposed.
c (ii) This update is postponed after the decomposition
c because of numerical accuracy and efficiency.
c
c the solution of the linear system is administered by subroutine solsy
c with calls to the Sylvester-Lyapunov solver we wrote.
c
c
c in addition to variables described previously, communication
c with prepj uses the following..
c y     = array containing predicted values on entry.
c wm    = real work space for matrices.  on output it contains the
c         decomposed blocks of the jacobian (kept in compact matrix
c         notation) and the accumulated orthogonal transformations
c         which we have used for the solution of the sylvester
c         equation and also eigenvalue information which has emerged
c         during the decomposition.
c         storage of matrix elements starts at wm(3).
c iwm   = integer work space containing eigenvalue information
c         starting at iwm(21). iwm also contains some stability
c         parameters ml = iwm(1) and mu = iwm(2).
c el0   = el(1) (input).
c ierpj = output error flag,  = 0 if no trouble, .gt. 0 if
c         p matrix found to be singular.
c jcur  = output flag = 1 to indicate that the jacobian matrix
c         is now current.
c this routine also uses the common variables el0, h, tn, uround,
c miter, n, nfe, and nje, and several from the block /intdre/.
c-----------------------------------------------------------------------
c next 2 lines added for the superstability check
      tcj = tn
      if (jstart .eq. 0) tpj = tcj
      nje = nje + 1
      ierpj = 0
      jcur = 1
      hl0 = h*el0
      lenp=2*(nqs+nps)+3*nn
      do 110 i = 1,lenp
 110    wm(i+2) = 0.0e0
      call probl(tn,y(neqn1),rarr, iarr)
      call doblks(nn,nnq,nnp,y(neqn1),y(neqn11),y(neqn12),
     *            y(neqn21),y(neqn22))
      call jac(nnq,nnp,y,y(neqn11),y(neqn12),y(neqn22),icontr)
      con = -hl0
c we have the "needed" blocks of Sylvester Equation in y(neqn11),y(neqn22)
c now, we first decompose them with calls to the Sylvester (Lyapunov)
c solver, then after we can multiply them by con and get the right
c signs.  Then, add identity (1/2 on each block)
      do 115 i=1,nps
 115     wm(2+i)=y(neqn22+i-1)
      do 120 i=1,nqs
 120     wm(2+2*nps+i)=y(neqn11+i-1)
c here we perform the orthogonal decomposition of the appropriate
c blocks of the Sylvester equation.  We differentiate depending on
c the possibility of symmetric DRE.
      if (icontr.eq.0) goto 275
      call axpxbd(wm(3),wm(np3),nnp,nnp,nnp,wm(np23),wm(np23q),nnq,
     *            nnq,nnq,uround,uround,wm(ntt),wm(nttp),wm(ntt2p),
     *            wm(ntt1),wm(ntt1q),wm(ntt2q),iwm(21),iwm(nnp+21),ier)
      if (ier.eq.4) then
         ierpj=1
         return
      endif
      do 125 i=1,nps
 125     wm(2+i)=wm(2+i)*con
      do 130 i=1,nqs
 130     wm(2+2*nps+i)=-wm(2+2*nps+i)*con
      j=3
      npp1=nnp+1
      do 135 i=1,nnp
         wm(j)=wm(j)+.5e0
 135     j=j+npp1
      j=3+2*nps
      nnqp1=nnq+1
      do 140 i=1,nnq
         wm(j)=wm(j)+.5e0
 140     j=j+nnqp1
      return
 275  continue
c load the correct matrix to be decomposed for the Lyapunov equation
      do 277 i=1,nps
 277     wm(2+i)=wm(2+2*nps+i)
      call atxxad(wm(3),wm(np3),nnp,nnp,nnp,uround,wm(ntt),wm(nttp),
     *            wm(ntt2p),iwm(21),ier)
      if (ier.eq.4) then
         ierpj=1
         return
      endif
c now complete the transformation by multiplying by con, getting the
c right sign and adding the identity (1/2, for us)
      do 280 i=1,nps
 280     wm(2+i)=-wm(2+i)*con
      j=3
      npp1=nnp+1
      do 285 i=1,nnp
         wm(j)=wm(j)+.5e0
 285     j=j+npp1
      return
c----------------------- end of subroutine prepj -----------------------
      end
      subroutine solsy (wm, x)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
c history: a routine by the same name is in the original lsode of
c          Alan Hindmarsh (LLL), but this is basically all new
c
      real wm, x
      dimension wm(*), x(*)
      integer iownd, iowns,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      integer nnq, nnp, icontr, nonlin, nsuper, nn, nns, nqs, nps,
     1   neqn1, neqn11, neqn12, neqn21, neqn22, np3, np23, np23q,
     2   ntt, nttp, ntt2p, ntt1, ntt1q, ntt2q, nflevc, nflevp
      real rowns,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      real reuvct, reuvpt, tcj, tpj
      common /ls0001/ rowns(209),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     3   iownd(14), iowns(6),
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      common /intdre/ reuvct, reuvpt, tcj, tpj, nnq, nnp, icontr,
     1   nonlin, nsuper, nn, nns, nqs, nps, neqn1, neqn11, neqn12,
     2   neqn21, neqn22, np3, np23, np23q, ntt, nttp, ntt2p, ntt1,
     3   ntt1q, ntt2q, nflevc, nflevp
c-----------------------------------------------------------------------
c this routine manages the solution of the linear system arising from
c a chord or quasi-Newton iteration.  it is called with miter=1.
c it calls the sylvester-lyapunov equation solver to do this
c
c communication with solsy uses the following variables..
c
c wm    = real work space containing the inverse diagonal matrix if
c         miter = 3 and the lu decomposition of the matrix otherwise.
c         storage of matrix elements starts at wm(3).
c x     = the right-hand side vector on input, and the solution vector
c         on output, of length n.
c iersl = output flag (in common).  iersl = 0 if no trouble occurred.
c *** we set iersl=-1 if troubles arose when solving Sylvester equation ***
c this routine also uses some of the common variables of /intdre/.
c-----------------------------------------------------------------------
      iersl = 0
      if (icontr.eq.0) goto 150
      call axpxbs(wm(3),wm(np3),nnp,nnp,nnp,wm(np23),wm(np23q),nnq,
     *            nnq,nnq,x,nnp,wm(ntt),wm(ntt1),iersl)
      if (iersl.eq.4) iersl=-1
      return
 150  call atxxas(wm(3),wm(np3),x,nnp,nnp,nnp,nnp,wm(ntt),iersl)
      if (iersl.eq.4) iersl=-1
      return
c----------------------- end of subroutine solsy -----------------------
      end
c
      subroutine jac(nnq,nnp,x,a011,a012,a022,icontr)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
      integer nnq, nnp, icontr
      real x, a011, a012, a022
      dimension a011(nnq,nnq),a012(nnq,nnp),a022(nnp,nnp),x(nnp,nnq)
c
c here we form the blocks of the Newton iteration (Sylvester eq.n)
c
      integer i, j, k
      real sum
      do 10 j=1,nnq
         do 10 i=1,nnq
              sum=0.e0
            do 5 k=1,nnp
  5            sum=sum+a012(i,k)*x(k,j)
 10         a011(i,j)=a011(i,j)+sum
      if (icontr.eq.0) goto 25
      do 20 j=1,nnp
          do 20 i=1,nnp
            sum=0.e0
            do 15 k=1,nnq
 15            sum=sum+x(i,k)*a012(k,j)
 20         a022(i,j)=a022(i,j)-sum
      return
 25   do 30 j=1,nnp
           do 30 i=1,nnp
 30           a022(i,j)=-a011(j,i)
      return
      end
c
      subroutine axpxbd(a,u,m,na,nu,b,v,n,nb,nv,epsa,
     1                 epsb,orta,era,eia,ortb,erb,eib,itypea,
     2                 itypeb,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart, Algorithm
c          432, Comm ACM 15, 1972.
c
      integer m, na, nu, n, nb, nv, itypea, itypeb, iflag
      real a, u, b, v, epsa, epsb, orta, era, eia, ortb, erb, eib
      dimension itypea(m), itypeb(n),
     1   a(na,m),u(nu,m),b(nb,n),v(nv,n),orta(m),ortb(n),era(m),
     2   eia(m),erb(n),eib(n)
c
c axpxbd is a fortran 77 subroutine to perform the decomposition of a and b
c for the equation ax + xb = c.  the matrices a and b are transformed
c into real schur form.
c the program requires the auxiliary
c subroutines orthes, ortran, hqr3.
c the parameters in the argument list are
c
c         a         a doubly subscripted array containing
c                   the matrix a. On return, the lower triangle
c                   and superdiagonal of the array a contain
c                   a lower real schur form of a.  the array
c                   a must be dimensioned at least m by m.
c         u         a doubly subscripted array that, on
c                   return, contains the orthogonal matrix
c                   that reduces a to real schur form with
c                   eigenvalues ordered from smallest to largest
c                   real parts.
c         m         the order of the matrix a.
c         na        the first dimension of the array a.
c         nu        the first dimension of the array u.
c         b         a doubly subscripted array containing
c                   the matrix b.  on return, the upper triangle
c                   and subdiagonal of the array b contain an
c                   upper real schur form of b.  the array b
c                   must be dimensioned at least n by n.
c         v         a doubly subscripted array that, on
c                   return, contains the orthogonal matrix
c                   that reduces b to real schur form with
c                   eigenvalues ordered from smallest to largest
c                   real parts.
c         n         the order of the matrix b.
c         nb        the first dimension of the array b.
c         nv        the first dimension of the array v.
c         epsa      a convergence criterion for the reduction
c                   of a to schur form.  epsa should be set
c                   to the unit roundoff if a is inputted
c                   with as many significant digits as possible.
c         epsb      a convergence criterrion for the reduction
c                   of b to real schur form.
c         iflag     an error indicator that, on return,
c                   contains an error signal.  if iflag is
c                   equal to 4 then the program was
c                   unable to reduce a (b) to real schur form.
c                   otherwise, iflag is untouched.
c         orta      a singly subscripted real array of dimension m
c                   needed as work space.
c         era       a singly subscripted real array of dimension m
c                   that on return contains the real parts of the
c                   eigenvalues of a.
c         eia       a singly subscripted real array of dimension m
c                   that on return contains the imaginary parts of
c                   the eigenvalues of a.
c         itypea    a singly subscripted integer array of dimension m
c                   whose i-th entry on return is
c                   0  i-th a-e.value is real
c                   1  i-th a-e.value is complex with positive imaginary part
c                   2  i-th a-e.value is complex with negative imaginary part
c                  -1  i-th a-e.value not successfully computed
c         ortb      same as orta, but for matrix b.
c         erb       same as era, but for the matrix b
c         eib       same as eia, but for matrix b
c         itypeb    same as itypea, but for the matrix b
c
c
c when epsa is negative the reduction of a to real schur form
c is skipped and the arrays a and u are assumed to contain the
c schur form and accompanying orthogonal matrix. likewise, if epsb
c is negative, the reduction of b to real schur form is skipped.
c this is ignored for our purposes of dresol.
c
      integer i, j, nlow, nup, mlow, mup
      real temp
c
      iflag=0
      mlow=1
      mup=m
      nlow=1
      nup=n
c
c if required,  reduce a(transpose) to upper real schur form.
c
      if (epsa .lt. 0.) go to 35
      do 10 i=1,m
         do 10 j=i,m
            temp = a(i,j)
            a(i,j) = a(j,i)
            a(j,i) = temp
   10 continue
      call orthes(na,m,mlow,mup,a,orta)
      call ortran(na,m,mlow,mup,a,orta,u)
      call hqr3(a,u,m,mlow,mup,epsa,era,eia,itypea,na,nu)
      do 20 i=mlow,mup
         if (itypea(i).lt.0) then
            iflag=4
            return
         endif
 20   continue
      do 30 i=1,m
         do 30 j=i,m
            temp = a(i,j)
            a(i,j) = a(j,i)
            a(j,i) = temp
   30 continue
c
c if required,  reduce b to upper real schur form.
c
   35 if (epsb .lt. 0.) goto 45
      call orthes(nb,n,nlow,nup,b,ortb)
      call ortran(nb,n,nlow,nup,b,ortb,v)
      call hqr3(b,v,n,nlow,nup,epsb,erb,eib,itypeb,nb,nv)
      do 40 i=nlow,nup
         if (itypeb(i).lt.0) then
            iflag=4
            return
         endif
 40   continue
 45   return
c end of axpxbd
      end
c
      subroutine axpxbs(a,u,m,na,nu,b,v,n,nb,nv,c,nc,orta,ortb,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart.
c
c axpxbs is a fortran 77 subroutine to solve the sylvester
c equation ax + xb = c.  the matrices a and b are already in
c real schur form.
c the program requires the auxiliary
c subroutines shrslv and sysslv.
c the parameters in the argument list are
c
c         a         a doubly subscripted array as follow.
c                   the lower triangle
c                   and superdiagonal of the array a contain
c                   a lower real schur form of a.  the array
c                   a must be dimensioned at least m by m.
c         u         a doubly subscripted array that,
c                   contains the orthogonal matrix
c                   that has reduced a to real schur form with
c                   eigenvalues ordered from smallest to
c                   largest real parts.
c         m         the order of the matrix a.
c         na        the first dimension of the array a.
c         nu        the first dimension of the array u.
c         b         a doubly subscripted array as follow.
c                   the upper triangle
c                   and subdiagonal of the array b contain an
c                   upper real schur form of b.  the array b
c                   must be dimensioned at least n by n.
c         v         a doubly subscripted array that
c                   contains the orthogonal matrix
c                   that has reduced b to real schur form with
c                   eigenvalues ordered from smallest to
c                   largest real parts.
c         n         the order of the matrix b.
c         nb        the first dimension of the array b.
c         nv        the first dimension of the array v.
c         c         a doubly subscripted array that contains
c                   the right-hand-side matrix c.  on return
c                   it contains the solution matrix x.
c         nc        the first dimension of c.
c         iflag     an error indicator that, on return,
c                   contains an error signal.  if iflag is
c                   equal to 4 then the program detected a division
c                   by zero.  this means sylvester equation not
c                   solvable.  otherwise, iflag is set to 0.
c         orta      a singly subscripted real array of dimension m
c                   needed as work space.
c         ortb      same as orta, but of dimension n.
c
c
      integer
     *na,nu,nb,nv,nc,m,n,i,j,k,iflag
      real
     *a(na,m),u(nu,m),b(nb,n),v(nv,n),orta(m),ortb(n),c(nc,n)
      iflag=0
c transform c.
      do 20 j=1,n
         do 10 i=1,m
            orta(i)=0.
            do 10 k=1,m
               orta(i) = orta(i) + u(k,i)*c(k,j)
   10 continue
      do 20 i=1,m
         c(i,j) = orta(i)
   20 continue
      do 40 i=1,m
         do 30 j=1,n
            ortb(j) = 0.
            do 30 k=1,n
               ortb(j) = ortb(j) + c(i,k)*v(k,j)
   30    continue
         do 40 j=1,n
            c(i,j) = ortb(j)
   40 continue
c
c solve the transformed system.
c
      call shrslv(a,b,c,m,n,na,nb,nc,iflag)
      if (iflag.eq.4) return
c
c transform c back to the solution.
c
      do 100 j=1,n
         do 90 i=1,m
            orta(i) = 0.
            do 90 k=1,m
               orta(i) = orta(i) + u(i,k)*c(k,j)
   90 continue
      do 100 i=1,m
         c(i,j) = orta(i)
  100 continue
      do 120 i=1,m
         do 110 j=1,n
            ortb(j) = 0.
            do 110 k=1,n
               ortb(j) = ortb(j) + c(i,k)*v(j,k)
  110 continue
      do 120 j=1,n
         c(i,j) = ortb(j)
  120 continue
      return
      end
c
      subroutine shrslv(a,b,c,m,n,na,nb,nc,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart.
c
c shrslv is a fortran iv subroutine to solve the real matrix
c equation ax + bx = c,where a is in lower real schur form
c and b is in upper real schur form.  shrslv uses the aux-
c iliary subroutine sysslv,which it communicates with
c through the common block slvblk.  the parameters in the
c argument list are
c           a      a doubly subscripted array containing the
c                  matrix a in lower real schur form.
c           b      a doubly subscripted array containing the
c                  matrix b in upper real schur form.
c           c      a doubly subscripted array containing rhe
c                  matrix c.
c           m      the order of the matrix a.
c           n      the order of the matrix b.
c           na     the first demension of the array a.
c           nb     the first dimension of the array b.
c           nc     the first dimension of the array c.
c
      integer
     1m,n,na,nb,nc,k,km1,dk,kk,l,lm1,dl,ll,i,ib,j,ja,nsys,iflag
      real
     1a(na,m),b(nb,n),c(nc,n),t,p
      common/slvblk/t(5,5),p(5),nsys
      l = 1
   10   lm1 = l-1
        dl = 1
        if (l .eq. n) go to 15
        if (b(l+1,l) .ne. 0.) dl = 2
   15   ll = l+dl-1
        if (l .eq. 1) go to 30
        do 20 j=l,ll
           do 20 i=1,m
              do 20 ib=1,lm1
                 c(i,j) = c(i,j) - c(i,ib)*b(ib,j)
   20   continue
   30   k = 1
   40     km1 = k-1
          dk = 1
          if (k .eq. m) go to 45
          if (a(k,k+1) .ne. 0.) dk = 2
   45     kk = k+dk-1
          if (k .eq. 1) go to 60
          do 50 i=k,kk
             do 50 j=l,ll
                do 50 ja=1,km1
                   c(i,j) = c(i,j) - a(i,ja)*c(ja,j)
   50     continue
   60     if (dl .eq. 2) go to 80
          if (dk .eq. 2) go to 70
          t(1,1) = a(k,k) + b(l,l)
          if (t(1,1) .eq. 0.) then
             iflag=4
             return
          endif
          c(k,l) = c(k,l)/t(1,1)
          go to 100
   70     t(1,1) = a(k,k) + b(l,l)
          t(1,2) = a(k,kk)
          t(2,1) = a(kk,k)
          t(2,2) = a(kk,kk) + b(l,l)
          p(1) = c(k,l)
          p(2) = c(kk,l)
          nsys = 2
          call sysslv
          c(k,l) = p(1)
          c(kk,l) = p(2)
          go to 100
   80     if (dk .eq. 2) go to 90
          t(1,1) = a(k,k) + b(l,l)
          t(1,2) = b(ll,l)
          t(2,1) = b(l,ll)
          t(2,2) = a(k,k) + b(ll,ll)
          p(1) = c(k,l)
          p(2) = c(k,ll)
          nsys = 2
          call sysslv
          c(k,l) = p(1)
          c(k,ll) = p(2)
          go to 100
   90     t(1,1) = a(k,k) + b(l,l)
          t(1,2) = a(k,kk)
          t(1,3) = b(ll,l)
          t(1,4) = 0.
          t(2,1) = a(kk,k)
          t(2,2) = a(kk,kk) + b(l,l)
          t(2,3) = 0.
          t(2,4) = t(1,3)
          t(3,1) = b(l,ll)
          t(3,2) = 0.
          t(3,3) = a(k,k) + b(ll,ll)
          t(3,4) = t(1,2)
          t(4,1) = 0.
          t(4,2) = t(3,1)
          t(4,3) = t(2,1)
          t(4,4) = a(kk,kk) + b(ll,ll)
          p(1) = c(k,l)
          p(2) = c(kk,l)
          p(3) = c(k,ll)
          p(4) = c(kk,ll)
          nsys = 4
          call sysslv
          c(k,l) = p(1)
          c(kk,l) = p(2)
          c(k,ll) = p(3)
          c(kk,ll) = p(4)
  100   k = k + dk
        if (k .le. m) go to 40
      l = l + dl
      if (l .le. n) go to 10
      return
      end
c
c
      subroutine atxxad(a,u,n,na,nu,eps,ort,er,ei,itype,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart.
c
c
c atxxad is a fortran 77 subroutine to decompose a prior to solving the
c real lyapunov equation trans(a)*x + x*a = c, where c is symmetric and
c trans(a) denotes the transpose of a. the equation is
c transformed so that a is in upper real schur form.
c the program requires the auxiliary subroutines orthes,
c ortan, hqr3 and related. the parameters in the argument
c list are
c           a     a doubly subscripted array containing the
c                 matrix a. on return, the upper triangle
c                 and the first subdiagonal of the array a
c                 contain an upper real schur form of a.
c                 the array a must be dimensioned ar least
c                 n by n.
c           u     a doubly subscripted array that, on
c                 return, contains the orthogonal matrix
c                 that reduces a to upper real schur form with
c                 eigenvalues ordered in descending absolute value.
c           n     the order of the matrix a.
c           na    the first dimension of the the array a.
c           nu    the first dimension of the array u.
c           eps   a convergence criterion for the reduction
c                 of a to real schur form. eps should be
c                 set to the machine precision.
c                 when a is correct within the precision itself.
c           ort   singly subscripted real array of dimension n
c                 needed as work space
c           er    singly subscripted real array of dimension n
c                 that on return contains the real parts of the e.values
c           ei    singly subscripted real array that on return
c                 contains the imaginary parts of the e.values
c         itype   integer array whose i-th entry on return is
c                 0  i-th e.value is real
c                 1  i-th e.value is complex with positive imaginary part
c                 2  i-th e.value is complex with negative imaginary part
c                -1  i-th e.value not successfully computed
c          iflag  an error indicator.  if iflag is set to 4, then
c                 program was unable to reduce to a to real schur
c                 form. it is set to 0 otherwise.
c
c if eps is negative, the reduction of a to real schur
c form is skipped and the arrays a and u are assumed to
c contain the schur form and accompanying orthogonal matrix.
c this makes sense just to check the input.
c
      integer
     *i,n,na,nu,itype(n),nlow,nup,iflag
      real
     *a(na,n),u(nu,n),eps,ort(n),er(n),ei(n)
      iflag=0
      nlow=1
      nup=n
c
c if required,  reduce a to upper real schur form.
c
      if (eps .lt. 0.) go to 15
      call orthes(na,n,nlow,nup,a,ort)
      call ortran(na,n,nlow,nup,a,ort,u)
      call hqr3(a,u,n,nlow,nup,eps,er,ei,itype,na,nu)
      do 10 i=nlow,nup
         if (itype(i).lt.0)then
            iflag=4
            return
         endif
 10   continue
 15   return
c end of atxxad
      end
c
      subroutine atxxas(a,u,c,n,na,nu,nc,ort,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart.
c
c atxxas is a fortran iv subroutine to solve the real matrix
c equation trans(a)*x + x*a = c, where c is symmetric and
c trans(a) denotes the transpose of a.
c a is in upper real schur form, and the
c transformed equation is solved by a recursive procedure.
c the program requires the auxiliary subroutines symslv and sysslv
c
c parameters in the calling list are>
c           a     a doubly subscripted array as follows.
c                 the upper triangle
c                 and the first subdiagonal of the array a
c                 contain an upper real schur form of a.
c                 the array a must be dimensioned ar least
c                 n by n.
c           u     a doubly subscripted array that
c                 contains the orthogonal matrix
c                 that has reduced a to upper real schur form with
c                 eigenvalues ordered in descending absolute value.
c           c     a doubly subscripted array containing
c                 the matrix c. on return, c contains
c                 the solution matrix x.
c           n     the order of the matrix a.
c           na    the first dimension of the the array a.
c           nu    the first dimension of the array u.
c           nc    the first dimension of the array c.
c           ort   singly subscripted real array of dimension n
c                 needed as work space
c          iflag  an error indicator.  if iflag is set to 4, then
c                 program detected a zero eigenvalue of a, thus
c                 cannot solve the lyapunov equation. set to 0 otherwise.
c
      integer
     *i,j,k,n,na,nu,nc,iflag
      real
     *a(na,n),u(nu,n),c(nc,n),ort(n)
      iflag=0
c transform c.
      do 20 i=1,n
         c(i,i) = c(i,i)/2.
   20 continue
      do 40 i=1,n
         do 30 j=1,n
            ort(j)=0.
            do 30 k=i,n
               ort(j) = ort(j) + c(i,k)*u(k,j)
   30 continue
      do 40 j=1,n
         c(i,j) = ort(j)
   40 continue
      do 60 j=1,n
         do 50 i=1,n
            ort(i) = 0.
            do 50 k=1,n
               ort(i) = ort(i) + u(k,i)*c(k,j)
   50 continue
      do 60 i=1,n
         c(i,j) = ort(i)
   60 continue
      do 70 i=1,n
         do 70 j=i,n
            c(i,j)=c(i,j)+c(j,i)
            c(j,i)=c(i,j)
   70 continue
c
c solve the transformed system.
c
      call symslv(a,c,n,na,nc,iflag)
      if (iflag.eq.4) return
c
c transform c back to the solution.
c
      do 80 i=1,n
         c(i,i)=c(i,i)/2.
   80 continue
         do 100 i=1,n
            do 90 j=1,n
            ort(j) = 0.
            do 90 k=i,n
               ort(j) = ort(j) + c(i,k)*u(j,k)
   90 continue
      do 100 j=1,n
         c(i,j) = ort(j)
  100 continue
      do 120 j=1,n
         do 110 i=1,n
            ort(i) = 0.
            do 110 k=1,n
               ort(i) = ort(i) + u(i,k)*c(k,j)
  110 continue
      do 120 i=1,n
         c(i,j) = ort(i)
  120 continue
      do 130 i=1,n
         do 130 j=i,n
            c(i,j)=c(i,j)+c(j,i)
            c(j,i)=c(i,j)
  130 continue
      return
      end
c
c
      subroutine symslv(a,c,n,na,nc,iflag)
c
c last modification:  august 12, 1989
c
c history: based on the software by Bartels & Stewart.
c
c
c symslv is fortran 77 subroutine to solve the real matrix
c equation trans(a)*x + xa = c, where c is symmetric, a is
c in upper real schur form, and trans(a) denotes the transpose
c of a. symslv uses the auxiliary subroutine sysslv, which it
c communicates with through the common block slvblk. the parameters
c in the argument list are
c        a      a doubly subscripted array containing the matrix a
c               in upper real schur form.
c        c      a doubly subscripted array containing the matrix c.
c        n      the order of the matrix a.
c        na     the first dimension of the array a.
c        nc     the first dimension of the array c.
c
      integer
     *n,na,nc,k,km1,dk,kk,l,ldl,dl,ll,i,ia,j,nsys,iflag
      real
     *a(na,n),c(nc,n),t,p
      common/slvblk/t(5,5),p(5),nsys
      l = 1
   10   dl = 1
        if (l .eq. n) go to 20
        if (a(l+1,l) .ne. 0.) dl = 2
   20   ll = l+dl-1
        k = l
   30     km1 = k-1
          dk = 1
          if (k .eq. n) go to 35
          if (a(k+1,k) .ne. 0.) dk = 2
   35     kk = k+dk-1
          if (k .eq. l) go to 45
          do 40 i=k,kk
             do 40 j=l,ll
                do 40 ia=l,km1
                   c(i,j) = c(i,j) - a(ia,i)*c(ia,j)
   40     continue
   45     if (dl .eq. 2) go to 60
          if (dk .eq. 2) go to 50
          t(1,1) = a(k,k) + a(l,l)
          if (t(1,1) .eq. 0.) then
             iflag=4
             return
          endif
          c(k,l) = c(k,l)/t(1,1)
          go to 90
   50     t(1,1) = a(k,k) + a(l,l)
          t(1,2) = a(kk,k)
          t(2,1) = a(k,kk)
          t(2,2) = a(kk,kk) + a(l,l)
          p(1) = c(k,l)
          p(2) = c(kk,l)
          nsys = 2
          call sysslv
          c(k,l) = p(1)
          c(kk,l) = p(2)
          go to 90
   60     if (dk .eq. 2) go to 70
          t(1,1) = a(k,k) + a(l,l)
          t(1,2) = a(ll,l)
          t(2,1) = a(l,ll)
          t(2,2) = a(k,k) + a(ll,ll)
          p(1) = c(k,l)
          p(2) = c(k,ll)
          nsys = 2
          call sysslv
          c(k,l) = p(1)
          c(k,ll) = p(2)
          go to 90
   70     if (k .ne. l) go to 80
          t(1,1) = a(l,l)
          t(1,2) = a(ll,l)
          t(1,3) = 0.
          t(2,1) = a(l,ll)
          t(2,2) = a(l,l)+a(ll,ll)
          t(2,3) = t(1,2)
          t(3,1) = 0.
          t(3,2) = t(2,1)
          t(3,3) = a(ll,ll)
          p(1) = c(l,l)/2.
          p(2) = c(ll,l)
          p(3) = c(ll,ll)/2.
          nsys = 3
          call sysslv
          c(l,l) = p(1)
          c(ll,l) = p(2)
          c(l,ll) = p(2)
          c(ll,ll) = p(3)
          go to 90
   80     t(1,1) = a(k,k) + a(l,l)
          t(1,2) = a(kk,k)
          t(1,3) = a(ll,l)
          t(1,4) = 0.
          t(2,1) = a(k,kk)
          t(2,2) = a(kk,kk) + a(l,l)
          t(2,3) = 0.
          t(2,4) = t(1,3)
          t(3,1) = a(l,ll)
          t(3,2) = 0.
          t(3,3) = a(k,k) + a(ll,ll)
          t(3,4) = t(1,2)
          t(4,1) = 0.
          t(4,2) = t(3,1)
          t(4,3) = t(2,1)
          t(4,4) = a(kk,kk) + a(ll,ll)
          p(1) = c(k,l)
          p(2) = c(kk,l)
          p(3) = c(k,ll)
          p(4) = c(kk,ll)
          nsys = 4
          call sysslv
          c(k,l) = p(1)
          c(kk,l) = p(2)
          c(k,ll) = p(3)
          c(kk,ll) = p(4)
   90   k = k + dk
        if (k .le. n) go to 30
        ldl = l + dl
        if (ldl .gt. n) return
        do 120 j=ldl,n
           do 100 i=l,ll
              c(i,j)=c(j,i)
  100      continue
           do 120 i=j,n
              do 110 k=l,ll
                 c(i,j)=c(i,j) - c(i,k)*a(k,j) - a(k,i)*c(k,j)
  110         continue
              c(j,i)=c(i,j)
  120      continue
      l=ldl
      go to 10
      end
c
      subroutine sysslv
c
c last modification:  august 12, 1989
c
c history: taken from the software by Bartels & Stewart.
c
c sysslv is a fortran 77 subroutine that solves the linear
c system ax = b of order n less than 5 by crout reduction
c followed by back substitution.  the matrix a,  the vector
c b,  and the order n are contained on the arrays a, b, and
c the variable n of the common block slvblk.  the solution
c is returned on the array b.
c
      integer
     *n,nm1,n1,k,km1,i,j,kp1,intr,im1,i1,ii
      real
     *rmax,a,b,aa,temp
      common/slvblk/a(5,5),b(5),n
      nm1 = n-1
      n1 = n+1
c
c compute the lu factorization of a.
c
      do 80 k=1,n
         km1 = k-1
         if (k .eq. 1) go to 20
         do 10 i=k,n
            do 10 j=1,km1
               a(i,k) = a(i,k) - a(i,j)*a(j,k)
   10    continue
   20    if (k .eq. n) go to 100
         kp1 = k+1
         rmax =abs(a(k,k))
         intr = k
         do 30 i=kp1,n
            aa = abs(a(i,k))
            if (aa .le. rmax) go to 30
            rmax = aa
            intr = i
   30    continue
         if (rmax .eq. 0.) stop
         a(n1,k) = intr
         if (intr .eq. k) go to 50
         do 40 j=1,n
            temp = a(k,j)
            a(k,j)  = a(intr,j)
            a(intr,j) = temp
   40    continue
   50    do 80 j=kp1,n
            if (k .eq. 1) go to 70
            do 60 i=1,km1
               a(k,j) = a(k,j) - a(k,i)*a(i,j)
   60       continue
   70       a(k,j)  = a(k,j)/a(k,k)
   80 continue
c
c interchange the components of b.
c
  100 do 110 j=1,nm1
         intr = a(n1,j)
         if (intr .eq. j) go to 110
         temp = b(j)
         b(j) = b(intr)
         b(intr) = temp
  110 continue
c
c solve lx = b.
c
  200 b(1) = b(1)/a(1,1)
      do 220 i=2,n
         im1 = i-1
         do 210 j=1,im1
            b(i) = b(i) - a(i,j)*b(j)
  210    continue
         b(i) = b(i)/a(i,i)
  220 continue
c
c solve ux = b.
c
  300 do 310 ii=1,nm1
         i = nm1-ii+1
         i1 = i+1
         do 310 j=i1,n
            b(i) = b(i) - a(i,j)*b(j)
  310 continue
      return
      end
c
c
      subroutine orthes(nm,n,low,igh,a,ort)
c
c this is from the EISPACK collection
c
      integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
      real a(nm,n),ort(igh)
      real f,g,h,scale
c *** functions
      real sqrt,abs,sign
c
c     this subroutine is a translation of the algol procedure orthes,
c     num. math. 12, 349-368(1968) by martin and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
c
c     given a real general matrix, this subroutine
c     reduces a submatrix situated in rows and columns
c     low through igh to upper hessenberg form by
c     orthogonal similarity transformations.
c
c     on input-
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement,
c
c        n is the order of the matrix,
c
c        low and igh are integers determined by the balancing
c          subroutine  balanc.  if  balanc  has not been used,
c          set low=1, igh=n,
c
c        a contains the input matrix.
c
c     on output-
c
c        a contains the hessenberg matrix.  information about
c          the orthogonal transformations used in the reduction
c          is stored in the remaining triangle under the
c          hessenberg matrix,
c
c        ort contains further information about the transformations.
c          only elements low through igh are used.
c
c     questions and comments should be directed to B. S. Garbow,
c     applied mathematics division, argonne national laboratory
c
c     ------------------------------------------------------------------
c
      la = igh - 1
      kp1 = low + 1
      if (la .lt. kp1) go to 200
c
      do 180 m = kp1, la
         h = 0.0
         ort(m) = 0.0
         scale = 0.0
c     ********** scale column (algol tol then not needed) **********
         do 90 i = m, igh
   90    scale = scale + abs(a(i,m-1))
c
         if (scale .eq. 0.0) go to 180
         mp = m + igh
c     ********** for i=igh step -1 until m do -- **********
         do 100 ii = m, igh
            i = mp - ii
            ort(i) = a(i,m-1) / scale
            h = h + ort(i) * ort(i)
  100    continue
c
         g = -sign(sqrt(h),ort(m))
         h = h - ort(m) * g
         ort(m) = ort(m) - g
c     ********** form (i-(u*ut)/h) * a **********
         do 130 j = m, n
            f = 0.0
c     ********** for i=igh step -1 until m do -- **********
            do 110 ii = m, igh
               i = mp - ii
               f = f + ort(i) * a(i,j)
  110       continue
c
            f = f / h
c
            do 120 i = m, igh
  120       a(i,j) = a(i,j) - f * ort(i)
c
  130    continue
c     ********** form (i-(u*ut)/h)*a*(i-(u*ut)/h) **********
         do 160 i = 1, igh
            f = 0.0
c     ********** for j=igh step -1 until m do -- **********
            do 140 jj = m, igh
               j = mp - jj
               f = f + ort(j) * a(i,j)
  140       continue
c
            f = f / h
c
            do 150 j = m, igh
  150       a(i,j) = a(i,j) - f * ort(j)
c
  160    continue
c
         ort(m) = scale * ort(m)
         a(m,m-1) = scale * g
  180 continue
c
  200 return
c     ********** last card of orthes **********
      end
c
      subroutine ortran(nm,n,low,igh,a,ort,z)
c
c this is from the EISPACK collection
c
      integer i,j,n,kl,mm,mp,nm,igh,low,mp1
      real a(nm,igh),ort(igh),z(nm,n)
      real g
c
c     this subroutine is a translation of the algol procedure ortrans,
c     num. math. 16, 181-204(1970) by peters and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
c
c     this subroutine accumulates the orthogonal similarity
c     transformations used in the reduction of a real general
c     matrix to upper hessenberg form by  orthes.
c
c     on input-
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement,
c
c        n is the order of the matrix,
c
c        low and igh are integers determined by the balancing
c          subroutine  balanc.  if  balanc  has not been used,
c          set low=1, igh=n,
c
c        a contains information about the orthogonal trans-
c          formations used in the reduction by  orthes
c          in its strict lower triangle,
c
c        ort contains further information about the trans-
c          formations used in the reduction by  orthes.
c          only elements low through igh are used.
c
c     on output-
c
c        z contains the transformation matrix produced in the
c          reduction by  orthes,
c
c        ort has been altered.
c
c     questions and comments should be directed to B. S. Garbow,
c     applied mathematics division, argonne national laboratory
c
c     ------------------------------------------------------------------
c
c     ********** initialize z to identity matrix **********
      do 80 i = 1, n
c
         do 60 j = 1, n
   60    z(i,j) = 0.0
c
         z(i,i) = 1.0
   80 continue
c
      kl = igh - low - 1
      if (kl .lt. 1) go to 200
c     ********** for mp=igh-1 step -1 until low+1 do -- **********
      do 140 mm = 1, kl
         mp = igh - mm
         if (a(mp,mp-1) .eq. 0.0) go to 140
         mp1 = mp + 1
c
         do 100 i = mp1, igh
  100    ort(i) = a(i,mp-1)
c
         do 130 j = mp, igh
            g = 0.0
c
            do 110 i = mp, igh
  110       g = g + ort(i) * z(i,j)
c     ********** divisor below is negative of h formed in orthes.
c                double division avoids possible underflow **********
            g = (g / ort(mp)) / a(mp,mp-1)
c
            do 120 i = mp, igh
  120       z(i,j) = z(i,j) + g * ort(i)
c
  130    continue
c
  140 continue
c
  200 return
c     ********** last card of ortran **********
      end
c
c below, the routines hqr3, exchng, qrstep and split
c are all from the HQR3 software of G.W.Stewart.
c minor modifications have been made to this software, i.e.
c translating it to single precision and little else
c
      subroutine hqr3 (a,v,n,nlow,nup,eps,er,ei,itype,na,nv)
c
c this routine is from the HQR3 software by G.W.Stewart
c
c modified on march 1988 by luca dieci
c so that it performs the reordering from
c smallest to largest real parts of the eigenvalues.
c
c     *****parameters:
      integer n,na,nlow,nup,nv,itype(n)
      real a(na,n),ei(n),er(n),eps,v(nv,n)
c
c     *****local variables:
      logical fail
      integer i,it,l,mu,nl,nu
      real e1,e2,p,q,r,s,t,w,x,y,z
c
c     *****functions:
      real abs
c
c     *****subroutines called:
c     exchng,qrstep,split
c
c     ------------------------------------------------------------------
c
c     *****purpose:
c     this subroutine reduces the upper hessenberg matrix a to quasi-
c     triangular form by unitary similarity transformations. the
c     eigenvalues of a, which are contained in the 1 x 1 and 2 x 2
c     diagonal blocks of the reduced matrix, are ordered in descending
c     order of magnitude along the diagonal. the transformations are
c     accumulated in the array v.
c
c     *****parameter description:
c     on input:
c        na,nv            row dimensions of the arrays containing a and
c                         v, respectively, as declared in the calling
c                         program dimension statement;
c
c        a                n x n array containing the upper hessenberg
c                         matrix to be reduced;
c
c        n                order of the matrices a and v;
c
c        nlow,nup         a(nlow,nlow-1) and a(nup,1+nup) are assumed
c                         to be zero, and only rows nlow through nup
c                         and columns nlow through nup are transformed,
c                         resulting in the calculation of eigenvalues
c                         nlow through nup;
c
c        eps              a convergence criterion used to determine when
c                         a subdiagonal element of a is negligible.
c                         specifically, a(i+1,i) is regarded as
c                         negligible if abs(a(i+1),i)).le.eps*
c                         (abs(a(i+1,i+1))). this means that the final
c                         matrix returned by the program will be exactly
c                         similar to a + e where e is of order
c                         eps*norm(a), for any reasonably balanced norm
c                         such as the row-sum norm;
c
c        itype            an integer vector of length n whose
c                         i-th entry is
c                         0  if the i-th eigenvalue is real,
c                         1  if the i-th eigenvalue is complex with
c                            positive imaginary part,
c                         2  if the i-th eigenvalue is complex with
c                            negative imaginary part,
c                         -1 if the i-th eigenvalue was not calculated
c                            successfully.
c
c        on output:
c
c        a                n x n array containing the reduced, quasi-
c                         triangular matrix;
c
c        v                n x n array containing the reducing
c                         transformations to be multiplied;
c
c        er,ei            real scratch vectors of length n which on
c                         return contain the real and imaginary parts,
c                         respectively, of the eigenvalues.
c
c     ------------------------------------------------------------------
c
      do 10 i=nlow,nup
         itype(i)=-1
10    continue
      t=0.0
      nu=nup
20    if (nu.lt.nlow) go to 240
      it=0
30    continue
      l=nu
40    continue
      if (l.eq.nlow) go to 50
      if (abs(a(l,l-1)).le.eps*(abs(a(l-1,l-1))+abs(a(l,l))))
     +    go to 50
      l=l-1
      go to 40
50    continue
      x=a(nu,nu)
      if (l.eq.nu) go to 160
      y=a(nu-1,nu-1)
      w=a(nu,nu-1)*a(nu-1,nu)
      if (l.eq.nu-1) go to 100
      if (it.eq.30) go to 240
      if (it.ne.10 .and. it.ne.20) go to 70
      t=t+x
      do 60 i=nlow,nu
         a(i,i)=a(i,i)-x
60    continue
      s=abs(a(nu,nu-1))+abs(a(nu-1,nu-2))
      x=0.75*s
      y=x
      w=-0.4375*s**2
70    continue
      it=it+1
      nl=nu-2
80    continue
      z=a(nl,nl)
      r=x-z
      s=y-z
      p=(r*s-w)/a(nl+1,nl)+a(nl,nl+1)
      q=a(nl+1,nl+1)-z-r-s
      r=a(nl+2,nl+1)
      s=abs(p)+abs(q)+abs(r)
      p=p/s
      q=q/s
      r=r/s
      if (nl.eq.l) go to 90
      if (abs(a(nl,nl-1))*(abs(q)+abs(r)).le.eps*abs(p)*
     +   (abs(a(nl-1,nl-1))+abs(z)+abs(a(nl+1,nl+1)))) go to 90
      nl=nl-1
      go to 80
90    continue
      call qrstep (a,v,p,q,r,nl,nu,n,na,nv)
      go to 30
100   if (nu.ne.nlow+1) a(nu-1,nu-2)=0.0
      a(nu,nu)=a(nu,nu)+t
      a(nu-1,nu-1)=a(nu-1,nu-1)+t
      itype(nu)=0
      itype(nu-1)=0
      mu=nu
110   continue
      nl=mu-1
      call split (a,v,n,nl,e1,e2,na,nv)
      if (a(mu,mu-1).eq.0.0) go to 170
      if (mu.eq.nup) go to 230
      if (mu.eq.nup-1) go to 130
      if (a(mu+2,mu+1).eq.0.0) go to 130
      if (a(mu-1,mu-1)+a(mu,mu).le.a(mu+1,mu+1)+
     +   a(mu+2,mu+2)) go to 230
      call exchng (a,v,n,nl,2,2,eps,fail,na,nv)
      if (.not.fail) go to 120
      itype(nl)=-1
      itype(nl+1)=-1
      itype(nl+2)=-1
      itype(nl+3)=-1
      go to 240
120   continue
      mu=mu+2
      go to 150
130   continue
      if (a(mu-1,mu-1)+a(mu,mu).le.
     +   2.*a(mu+1,mu+1)) go to 230
      call exchng (a,v,n,nl,2,1,eps,fail,na,nv)
      if (.not.fail) go to 140
      itype(nl)=-1
      itype(nl+1)=-1
      itype(nl+2)=-1
      go to 240
140   continue
      mu=mu+1
150   continue
      go to 110
160   nl=0
      a(nu,nu)=a(nu,nu)+t
      if (nu.ne.nlow) a(nu,nu-1)=0.0
      itype(nu)=0
      mu=nu
170   continue
180   continue
      if (mu.eq.nup) go to 220
      if (mu.eq.nup-1) go to 200
      if (a(mu+2,mu+1).eq.0.0) go to 200
      if (2.*a(mu,mu).le.a(mu+1,mu+1)+a(mu+2,mu+2))
     +    go to 220
      call exchng (a,v,n,mu,1,2,eps,fail,na,nv)
      if (.not.fail) go to 190
      itype(mu)=-1
      itype(mu+1)=-1
      itype(mu+2)=-1
      go to 240
190   continue
      mu=mu+2
      go to 210
200   continue
      if (a(mu,mu).le.a(mu+1,mu+1)) go to 220
      call exchng (a,v,n,mu,1,1,eps,fail,na,nv)
      mu=mu+1
210   continue
      go to 180
220   continue
      mu=nl
      nl=0
      if (mu.ne.0) go to 170
230   continue
      nu=l-1
      go to 20
240   if (nu.lt.nlow) go to 260
      do 250 i=nlow,nu
         a(i,i)=a(i,i)+t
250   continue
260   continue
      nu=nup
270   continue
      if (itype(nu).ne.(-1)) go to 280
      nu=nu-1
      go to 310
280   continue
      if (nu.eq.nlow) go to 290
      if (a(nu,nu-1).eq.0.0) go to 290
      call split (a,v,n,nu-1,e1,e2,na,nv)
      if (a(nu,nu-1).eq.0.0) go to 290
      er(nu)=e1
      ei(nu-1)=e2
      er(nu-1)=er(nu)
      ei(nu)=-ei(nu-1)
      itype(nu-1)=1
      itype(nu)=2
      nu=nu-2
      go to 300
290   continue
      er(nu)=a(nu,nu)
      ei(nu)=0.0
      nu=nu-1
300   continue
310   continue
      if (nu.ge.nlow) go to 270
      return
c
c     last line of hqr3
c
      end
      subroutine exchng (a,v,n,l,b1,b2,eps,fail,na,nv)
c
c
c     based on stewart's hqr3 software.
c     (hammarling's correction added by ajl 01/23/86)
c
c     *****parameters
      integer n,b1,b2,l,na,nv
      real a(na,n),eps,v(nv,n)
      logical fail
c
c     *****local variables:
      integer i,it,j,l1,m
      real p,q,r,s,w,x,y,z
c
c     *****functions:
      real abs,sqrt,amax1
c
c     *****subroutines called:
c     qrstep
c
c     ------------------------------------------------------------------
c
c     *****purpose:
c     given the upper schur matrix a with consecutive b1 x b1 and
c     b2 x b2 diagonal blocks (b1, b2.le.2) starting at a(l,l), this
c     subroutine produces a unitary similarity transformation that
c     exchanges the blocks along with their eigenvalues. the
c     transformation is accumulated in v.
c
c     *****parameter description:
c     on input:
c        na,nv            row dimensions of the arrays containing a
c                         and v, respectively, as declared in the
c                         calling program dimension statement;
c
c        a                n x n matrix whose blocks are to be
c                         interchanged;
c
c        n                order of the matrix a;
c
c        l                position of the blocks;
c
c        b1               an integer containing the size of the first
c                         block;
c
c        b2               an integer containing the size of the second
c                         block;
c
c        eps              a convergence criterion (cf. hqr3).
c
c     on output:
c
c        fail             a logical variable which is .false. on a
c                         normal return. if thirty iterations were
c                         performed without convergence, fail is set to
c                         .true. and the element a(l+b2,l+b2-1) cannot
c                         be assumed zero.
c
c     ------------------------------------------------------------------
c
      fail=.false.
      if (b1.eq.2) go to 70
      if (b2.eq.2) go to 40
      l1=l+1
      q=a(l+1,l+1)-a(l,l)
      if (q .eq. 0.0) return
      p=a(l,l+1)
      r=amax1(abs(p),abs(q))
      if (r.eq.0.0) r = q
      p=p/r
      q=q/r
      r=sqrt(p**2+q**2)
      p=p/r
      q=q/r
      do 10 j=l,n
         s=p*a(l,j)+q*a(l+1,j)
         a(l+1,j)=p*a(l+1,j)-q*a(l,j)
         a(l,j)=s
10    continue
      do 20 i=1,l1
         s=p*a(i,l)+q*a(i,l+1)
         a(i,l+1)=p*a(i,l+1)-q*a(i,l)
         a(i,l)=s
20    continue
      do 30 i=1,n
         s=p*v(i,l)+q*v(i,l+1)
         v(i,l+1)=p*v(i,l+1)-q*v(i,l)
         v(i,l)=s
30    continue
      a(l+1,l)=0.0
      return
40    continue
      x=a(l,l)
      p=1.0
      q=1.0
      r=1.0
      call qrstep (a,v,p,q,r,l,l+2,n,na,nv)
      it=0
50    it=it+1
      if (it.le.30) go to 60
      fail=.true.
      return
60    continue
      p=a(l,l)-x
      q=a(l+1,l)
      r=0.0
      call qrstep (a,v,p,q,r,l,l+2,n,na,nv)
      if (abs(a(l+2,l+1)).gt.eps*(abs(a(l+1,l+1))+abs(a(l+2,l+2))))
     +    go to 50
      a(l+2,l+1)=0.0
      return
70    continue
      m=l+2
      if (b2.eq.2) m=m+1
      x=a(l+1,l+1)
      y=a(l,l)
      w=a(l+1,l)*a(l,l+1)
      p=1.0
      q=1.0
      r=1.0
      call qrstep (a,v,p,q,r,l,m,n,na,nv)
      it=0
80    it=it+1
      if (it.le.30) go to 90
      fail=.true.
      return
90    continue
      z=a(l,l)
      r=x-z
      s=y-z
      p=(r*s-w)/a(l+1,l)+a(l,l+1)
      q=a(l+1,l+1)-z-r-s
      r=a(l+2,l+1)
      s=abs(p)+abs(q)+abs(r)
      p=p/s
      q=q/s
      r=r/s
      call qrstep (a,v,p,q,r,l,m,n,na,nv)
      if (abs(a(m-1,m-2)).gt.eps*(abs(a(m-1,m-1))+abs(a(m-2,m-2))))
     +    go to 80
      a(m-1,m-2)=0.0
      return
c
c     last line of exchng
c
      end
      subroutine qrstep (a,v,p,q,r,nl,nu,n,na,nv)
c
c     based on stewart's hqr3 software.
c
c     *****parameters:
      integer n,na,nl,nu,nv
      real a(na,n),p,q,r,v(nv,n)
c
c     *****local variables:
      logical last
      integer i,j,k,nl2,nl3,num1
      real s,x,y,z
c
c     *****functions:
      integer min0
      real abs,sqrt
c
c     *****subroutines called:
c     none
c
c     ------------------------------------------------------------------
c
c     *****purpose:
c     this subroutine performs one implicit qr step on the upper
c     hessenberg matrix a. the shift is determined by the numbers p,q,
c     and r, and the step is applied to rows and columns nl throuth nu.
c     the transformations are accumulated in the array v.
c
c     *****parameter description:
c     on input:
c        na,nv            row dimensions of the arrays containing a
c                         and v, respectively, as declared in the
c                         calling program dimension statement;
c
c        a                n x n upper hessenberg matrix on which the qr
c                         step is to be performed;
c
c        p,q,r            parameters which determine the shift;
c
c        nl               the lower limit of the step;
c
c        nu               the upper limit of the step;
c
c        n                order of the matrix a.
c
c     on output:
c
c        v                n x n real scratch array containing the
c                         accumulated transformations.
c     ------------------------------------------------------------------
c
      nl2=nl+2
      do 10 i=nl2,nu
         a(i,i-2)=0.0
10    continue
      if (nl2.eq.nu) go to 30
      nl3=nl+3
      do 20 i=nl3,nu
         a(i,i-3)=0.0
20    continue
30    continue
      num1=nu-1
      do 130 k=nl,num1
         last=k.eq.num1
         if (k.eq.nl) go to 40
         p=a(k,k-1)
         q=a(k+1,k-1)
         r=0.0
         if (.not.last) r=a(k+2,k-1)
         x=abs(p)+abs(q)+abs(r)
         if (x.eq.0.0) go to 130
         p=p/x
         q=q/x
         r=r/x
40       continue
         s=sqrt(p**2+q**2+r**2)
         if (p.lt.0.0) s=-s
         if (k.eq.nl) go to 50
         a(k,k-1)=-s*x
         go to 60
50       continue
         if (nl.ne.1) a(k,k-1)=-a(k,k-1)
60       continue
         p=p+s
         x=p/s
         y=q/s
         z=r/s
         q=q/p
         r=r/p
         do 80 j=k,n
            p=a(k,j)+q*a(k+1,j)
            if (last) go to 70
            p=p+r*a(k+2,j)
            a(k+2,j)=a(k+2,j)-p*z
70          continue
            a(k+1,j)=a(k+1,j)-p*y
            a(k,j)=a(k,j)-p*x
80       continue
         j=min0(k+3,nu)
         do 100 i=1,j
            p=x*a(i,k)+y*a(i,k+1)
            if (last) go to 90
            p=p+z*a(i,k+2)
            a(i,k+2)=a(i,k+2)-p*r
90          continue
            a(i,k+1)=a(i,k+1)-p*q
            a(i,k)=a(i,k)-p
100      continue
         do 120 i=1,n
            p=x*v(i,k)+y*v(i,k+1)
            if (last) go to 110
            p=p+z*v(i,k+2)
            v(i,k+2)=v(i,k+2)-p*r
110         continue
            v(i,k+1)=v(i,k+1)-p*q
            v(i,k)=v(i,k)-p
120      continue
130   continue
      return
c
c     last line of qrstep
c
      end
      subroutine split (a,v,n,l,e1,e2,na,nv)
c
c     based on stewart's hqr3 software.
c
c     *****parameters:
      integer l,n,na,nv
      real a(na,n),v(nv,n),e1,e2
c
c     *****local variables:
      integer i,j,l1
      real p,q,r,t,u,w,x,y,z
c
c     *****functions:
      real abs,sqrt
c
c     ****subroutines called:
c     none
c
c     ------------------------------------------------------------------
c
c     *****purpose:
c     given the upper-hessenberg matrix a with a 2 x 2 block starting at
c     a(l,l), this program determines if the corresponding eigenvalues
c     are real or complex. if they are real, a rotation is determined
c     that reduces the block to upper-triangular form with the
c     eigenvalue of largest absolute value appearing first. the
c     rotation is accumulated in the array v.
c
c     *****parameter description:
c     on input:
c        na,nv            row dimensions of the arrays containing
c                         a and v, respectively, as declared in the
c                         calling program dimension statement;
c
c        a                the upper hessenberg matrix whose 2 x 2 block
c                         is to be split;
c
c        n                order of the matrix a;
c
c        l                position of the 2 x 2 block.
c
c     on output:
c
c        v                an n x n array containing the accumulated
c                         splitting transformation;
c
c        e1,e2            real scalars. if the eigenvalues are complex,
c                         e1 and e2 contain their common real part and
c                         positive imaginary part (respectively).
c                         if the eigenvalues are real, e1 contains the
c                         one largest in absolute value and e2 contains
c                         the other one.
c
c     ------------------------------------------------------------------
c
      x=a(l+1,l+1)
      y=a(l,l)
      w=a(l,l+1)*a(l+1,l)
      p=(y-x)/2.0
      q=p**2+w
      if (q.ge.0.0) go to 10
      e1=p+x
      e2=sqrt(-q)
      return
10    continue
      z=sqrt(q)
      if (p.lt.0.0) go to 20
      z=p+z
      go to 30
20    continue
      z=p-z
30    continue
      if (z.eq.0.0) go to 40
      r=-w/z
      go to 50
40    continue
      r=0.0
50    continue
      if (abs(x+z).ge.abs(x+r)) z=r
      y=y-x-z
      x=-z
      t=a(l,l+1)
      u=a(l+1,l)
      if (abs(y)+abs(u).le.abs(t)+abs(x)) go to 60
      q=u
      p=y
      go to 70
60    continue
      q=x
      p=t
70    continue
      r=sqrt(p**2+q**2)
      if (r.gt.0.0) go to 80
      e1=a(l,l)
      e2=a(l+1,l+1)
      a(l+1,l)=0.0
      return
80    continue
      p=p/r
      q=q/r
      do 90 j=l,n
         z=a(l,j)
         a(l,j)=p*z+q*a(l+1,j)
         a(l+1,j)=p*a(l+1,j)-q*z
90    continue
      l1=l+1
      do 100 i=1,l1
         z=a(i,l)
         a(i,l)=p*z+q*a(i,l+1)
         a(i,l+1)=p*a(i,l+1)-q*z
100   continue
      do 110 i=1,n
         z=v(i,l)
         v(i,l)=p*z+q*v(i,l+1)
         v(i,l+1)=p*v(i,l+1)-q*z
110   continue
      a(l+1,l)=0.0
      e1=a(l,l)
      e2=a(l+1,l+1)
      return
c
c     last line of split
c
      end
c
      subroutine ewset (n, itol, rtol, atol, ycur, ewt)
c
c author: Alan Hindmarsh, LLL
c
c-----------------------------------------------------------------------
c this subroutine sets the error weight vector ewt according to
c     ewt(i) = rtol(i)*abs(ycur(i)) + atol(i),  i = 1,...,n,
c with the subscript on rtol and/or atol possibly replaced by 1 above,
c depending on the value of itol.
c-----------------------------------------------------------------------
      integer n, itol
      integer i
      real rtol, atol, ycur, ewt
      dimension rtol(*), atol(*), ycur(n), ewt(n)
c
      go to (10, 20, 30, 40), itol
 10   continue
      do 15 i = 1,n
 15     ewt(i) = rtol(1)*abs(ycur(i)) + atol(1)
      return
 20   continue
      do 25 i = 1,n
 25     ewt(i) = rtol(1)*abs(ycur(i)) + atol(i)
      return
 30   continue
      do 35 i = 1,n
 35     ewt(i) = rtol(i)*abs(ycur(i)) + atol(1)
      return
 40   continue
      do 45 i = 1,n
 45     ewt(i) = rtol(i)*abs(ycur(i)) + atol(i)
      return
c----------------------- end of subroutine ewset -----------------------
      end
      real function vnorm (n, v, w)
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
c history: this is a rewriting of function vnorm of lsode
c
c-----------------------------------------------------------------------
c this function routine computes the weighted max-norm
c of the vector of length n contained in the array v, with weights
c contained in the array w of length n..
c   vnorm = amax1( v(i)*w(i), i=1,n )
c-----------------------------------------------------------------------
      integer n, i
      real v, w, sum
c *** functions
      real amax1, abs
      dimension v(n), w(n)
      sum = 0.0e0
      do 10 i = 1,n
 10     sum = amax1(sum, abs(v(i)*w(i)))
      vnorm = sum
      return
c----------------------- end of function vnorm -------------------------
      end
      subroutine srcom (rsav, isav, job)
c
c revision: July 10 1990
c modifications' author: Luca Dieci, School of Mathematics, GaTech
c
c history: identical to the routine srcom of lsode by A.Hindmarsh,
c          but we added an extra block to save-restore
c
c-----------------------------------------------------------------------
c this routine saves or restores (depending on job) the contents of
c the common blocks ls0001, eh0001 and intdre which are used internally
c by one or more dresol solvers.
c
c rsav = real array of length 222 or more.
c isav = integer array of length 66 or more.
c job  = flag indicating to save or restore the common blocks..
c        job  = 1 if common is to be saved (written to rsav/isav)
c        job  = 2 if common is to be restored (read from rsav/isav)
c        a call with job = 2 presumes a prior call with job = 1.
c-----------------------------------------------------------------------
      integer isav, job
      integer ieh, ils, iint
      integer i, lenils, lenrls, lenrnt, lenint
      real rsav, rls, rint
      dimension rsav(*), isav(*)
      common /ls0001/ rls(218), ils(39)
      common /eh0001/ ieh(2)
      common /intdre/ rint(4), iint(25)
      data lenrls/218/, lenils/39/, lenrnt/4/, lenint/25/
c
      if (job .eq. 2) go to 100
c
      do 10 i = 1,lenrls
 10     rsav(i) = rls(i)
      do 15 i = 1,lenrnt
 15     rsav(lenrls+i) = rint(i)
      do 20 i = 1,lenils
 20     isav(i) = ils(i)
      isav(lenils+1) = ieh(1)
      isav(lenils+2) = ieh(2)
      do 25 i = 1,lenint
 25      isav(lenils+2+i) = iint(i)
      return
c
 100  continue
      do 110 i = 1,lenrls
 110     rls(i) = rsav(i)
      do 115 i = 1,lenrnt
 115     rint(i) = rsav(lenrls+i)
      do 120 i = 1,lenils
 120     ils(i) = isav(i)
      ieh(1) = isav(lenils+1)
      ieh(2) = isav(lenils+2)
      do 125 i = 1,lenint
 125     iint(i) = isav(lenils+2+i)
      return
c----------------------- end of subroutine srcom -----------------------
      end
      real function r1mach (idum)
      integer idum
c
c author: Alan Hindmarsh, LLL
c
c-----------------------------------------------------------------------
c this routine computes the unit roundoff of the machine.
c this is defined as the smallest positive machine number
c u such that  1.0 + u .ne. 1.0
c-----------------------------------------------------------------------
      real u, comp
      u = 1.0e0
 10   u = u*0.5e0
      comp = 1.0e0 + u
      if (comp .ne. 1.0e0) go to 10
      r1mach = u*2.0e0
      return
c----------------------- end of function r1mach ------------------------
      end
      subroutine xerrwv (msg, nmes, nerr, level, ni, i1, i2, nr, r1, r2)
      integer msg, nmes, nerr, level, ni, i1, i2, nr,
     1   i, lun, lunit, mesflg, ncpw, nch, nwds
      real r1, r2
      dimension msg(nmes)
c author: Alan Hindmarsh, LLL
c-----------------------------------------------------------------------
c subroutines xerrwv, xsetf, and xsetun, as given here, constitute
c a simplified version of the slatec error handling package.
c written by a. c. hindmarsh at llnl.  version of march 30, 1987.
c ----------------------------------------------------------------------
c
c all arguments are input arguments.
c
c msg    = the message (hollerith literal or integer array).
c nmes   = the length of msg (number of characters).
c nerr   = the error number (not used).
c level  = the error level..
c          0 or 1 means recoverable (control returns to caller).
c          2 means fatal (run is aborted--see note below).
c ni     = number of integers (0, 1, or 2) to be printed with message.
c i1,i2  = integers to be printed, depending on ni.
c nr     = number of reals (0, 1, or 2) to be printed with message.
c r1,r2  = reals to be printed, depending on nr.
c
c note..  this routine is machine-dependent and specialized for use
c in limited context, in the following ways..
c 1. the number of hollerith characters stored per word, denoted
c    by ncpw below, is a data-loaded constant.
c 2. the value of nmes is assumed to be at most 60.
c    (multi-line messages are generated by repeated calls.)
c 3. if level = 2, control passes to the statement   stop
c    to abort the run.  this statement may be machine-dependent.
c 4. r1 and r2 are assumed to be in single precision and are printed
c    in e21.13 format.
c-----------------------------------------------------------------------
c the following are instructions for installing this routine
c in different machine environments.
c
c to change the default output unit, change the data statement in the
c block data atlan2 below
c
c for a different number of characters per word, change the
c data statement setting ncpw below, and format 10.  alternatives for
c various computers are shown in comment cards.
c
c for a different run-abort command, change the statement following
c statement 100 at the end.
c-----------------------------------------------------------------------
      common /eh0001/ mesflg, lunit
c-----------------------------------------------------------------------
c the following data-loaded value of ncpw is valid for the cdc-6600
c and cdc-7600 computers.
c     data ncpw/10/
c the following is valid for the cray-1 and on the cyber 855 and 990.
c     data ncpw/8/
c the following is valid for the burroughs 6700 and 7800 computers.
c     data ncpw/6/
c the following is valid for the pdp-10 computer.
c     data ncpw/5/
c the following is valid for the vax computer with 4 bytes per integer,
c and for the ibm-360, ibm-370, ibm-303x, and ibm-43xx computers.
c and for the sun workstations sun3, sun4 and sparc
      data ncpw/4/
c the following is valid for the pdp-11, or vax with 2-byte integers.
c     data ncpw/2/
c-----------------------------------------------------------------------
c
      if (mesflg .eq. 0) go to 100
c get logical unit number. ---------------------------------------------
      lun = lunit
c get number of words in message. --------------------------------------
      nch = min0(nmes,60)
      nwds = nch/ncpw
      if (nch .ne. nwds*ncpw) nwds = nwds + 1
c write the message. ---------------------------------------------------
      write (lun, 10) (msg(i),i=1,nwds)
c-----------------------------------------------------------------------
c the following format statement is to have the form
c 10  format(1x,mmann)
c where nn = ncpw and mm is the smallest integer .ge. 60/ncpw.
c the following is valid for ncpw = 10.
c 10  format(1x,6a10)
c the following is valid for ncpw = 8.
c 10  format(1x,8a8)
c the following is valid for ncpw = 6.
c 10  format(1x,10a6)
c the following is valid for ncpw = 5.
c 10  format(1x,12a5)
c the following is valid for ncpw = 4.
  10  format(1x,15a4)
c the following is valid for ncpw = 2.
c 10  format(1x,30a2)
c-----------------------------------------------------------------------
      if (ni .eq. 1) write (lun, 20) i1
 20   format(6x,23hin above message,  i1 =,i10)
      if (ni .eq. 2) write (lun, 30) i1,i2
 30   format(6x,23hin above message,  i1 =,i10,3x,4hi2 =,i10)
      if (nr .eq. 1) write (lun, 40) r1
 40   format(6x,23hin above message,  r1 =,e21.13)
      if (nr .eq. 2) write (lun, 50) r1,r2
 50   format(6x,15hin above,  r1 =,e21.13,3x,4hr2 =,e21.13)
c abort the run if level = 2. ------------------------------------------
 100  if (level .ne. 2) return
      stop
c----------------------- end of subroutine xerrwv ----------------------
      end
      block data atlan2
c
c revision date: July 10 1990
c author: Luca Dieci, School of Mathematics, GaTech
c
c this block-data subprograms is needed to load the variables msflag and
c lunit into the common block eh0001.
c
c    the variables are as follows..
c       mesflg = print control flag..
c                1 means print all messages (the default).
c                0 means no printing.
c       lunit  = logical unit number for messages.
c                the default is 6 (machine-dependent).
      common /eh0001/ mesflg, lunit
      data mesflg/1/, lunit/6/
c----------------------- end of bloch data atlan2 ----------------------
      end
      subroutine xsetf (mflag)
c
c author: Alan Hindmarsh, LLL
c
c this routine resets the print control flag mflag.
c
      integer mflag, mesflg, lunit
      common /eh0001/ mesflg, lunit
c
      if (mflag .eq. 0 .or. mflag .eq. 1) mesflg = mflag
      return
c----------------------- end of subroutine xsetf -----------------------
      end
      subroutine xsetun (lun)
c
c author: Alan Hindmarsh, LLL
c
c this routine resets the logical unit number for messages.
c
      integer lun, mesflg, lunit
      common /eh0001/ mesflg, lunit
c
      if (lun .gt. 0) lunit = lun
      return
c----------------------- end of subroutine xsetun ----------------------
      end