      subroutine rowmap(n,f,ifcn,t,u,tend,hs,rtol,atol,itol,
     1                  jacv,ijacv,fdt,ifdt,solout,iout,work,
     2                  lwork,iwork,liwork,rpar,ipar,idid)
c-----!----------------------------------------------------------------- 
c
c   ROWMAP .. A ROW-code with Krylov techniques for large stiff ODEs.
c
c
c   ROWMAP solves the inital value problem for stiff systems of
c first order ODEs
c 
c      du     
c      -- = f(t,u),    u=(u(1),..,u(n)), f=(f(1),..f(n)).  
c      dt
c 
c
c   ROWMAP is based on the ROW-methods of order 4 of the code ROS4 of Hairer
c and Wanner (see [4]) and uses Krylov techniques for the solution of linear
c systems. By a special multiple Arnoldi process the order of the basic
c method is preserved with small Krylov dimensions.  Step size control is
c done by embedding with a method of order 3.
c
c Version of January 10, 1997.
c
c   Authors: 
c       H. Podhaisky, R. Weiner, 
c       Fachbereich Mathematik und Informatik, Universitaet Halle
c       06099 Halle, Germany
c       email: helmut@mail.mathematik.uni-halle.de
c              rw@mail.mathematik.uni-halle.de,
c     
c       B.A. Schmitt,
c       Fachbereich Mathematik, Universitaet Marburg
c       35032 Marburg, Germany  
c       email: schmitt@mathematik.uni-marburg.de
c
c
c-----!-----------------------------------------------------------------
c References:
c
c [1]  B.A. Schmitt and R. Weiner:
c         Matrix-free W-methods using a multiple Arnoldi Iteration,
c         APNUM 18(1995), 307-320
c
c [2]  R. Weiner, B.A. Schmitt an H. Podhaisky: 
c         ROWMAP - a ROW-code with Krylov techniques for large stiff
c         ODEs. Report 39, FB Mathematik und Informatik, 
c         Universitaet Halle, 1996
c       
c [3]  R.Weiner and B.A.Schmitt:
c         Consistency of Krylov-W-Methods in initial value problems,
c         Tech. Report 14, FB Mathematik und Informatik, 
c         Universitaet Halle, 1995
c
c [4]  E. Hairer and G. Wanner:
c         Solving Ordinary Differential Equations II, Springer-Verlag,
c         Second Edition, 1996
c   
c
c  For retrieving the latest version visit the ROWMAP World Wide Web 
c  homepage at URL :
c   http://www.mathematik.uni-halle.de/institute/numerik/software/
c
c-----!----------------------------------------------------------------- 
c
c   INPUT PARAMETERS
c   ----------------
c   
c   n      = Number of equations of the ODE system. 
c
c   f      = Name (external) of user-supplied subroutine computing the 
c            value f(t,u): 
c                subroutine f(n,t,u,udot,rpar,ipar)
c                real*8 t,u(n),udot(n),rpar(*)
c                integer n,ipar(*) 
c                udot(1)= ...
c            rpar,ipar (user parameters, see below).
c       
c   ifcn   = Switch:
c            ifcn = 0:  f(t,u) independent of t (autonomous)
c            ifcn = 1:  f(t,u) may depend on t (nonautonomous), default
c
c   t      = Initial t-value. Changed on exit.
c
c   u      = Initial values for u. Changed on exit.
c
c   tend   = Final t-value.
c
c   hs     = Initial step size guess. Changed on exit.
c            If hs=0d0 the code puts hs to a default value depending on 
c            atol, rtol and u. 
c           
c   rtol, atol = Relative and absolute error tolerances. They
c            can be both scalars or else both vectors of length n.
c            ROWMAP uses a weighted root-mean-square norm to measure the
c            size of error vectors. It is defined by
c
c                wrms=sqrt(sum[(err(i)/scal(i))**2, i=1,n]/n),
c
c            where
c                 scal(i)=atol+rtol*dabs(u(i))              (itol=0)
c            or
c                 scal(i)=atol(i)+rtol(i)*dabs(u(i))        (itol=1).
c
c   itol   = Switch for rtol,atol: 
c            itol=0: both are scalars (default).
c            itol=1: both are vectors.
c
c   jacv   = Name (external) of user-supplied subroutine computing 
c            the jacobian-vector product. 
c            z := Jac(t,u)*v, Jac=df/du:  
c                subroutine jacv(n,t,u,v,z,rpar,ipar)
c                real*8 t,u(n),v(n),z(n),rpar(*)
c                integer n,ipar(*)
c                z(1)=...   
c            rpar,ipar (see below).
c            This routine is only called if ijacv = 1. 
c            Supply a dummy subroutine in the case ijacv = 0
c            
c   ijacv  = Switch for the computation of the jacobian-vector products:
c            ijacv = 0: The products are approximated by finite differences.
c                       Subroutine jacv is never called (default).
c            ijacv = 1: The products are computed by the user-supplied 
c                       subroutine jacv.
c
c   fdt    = Name (external) of user-supplied subroutine computing the
c            partial derivate of f(t,u) with respect to t.
c            This routine is only called if idft=1 and ifcn=1.
c
c                subroutine fdt(n,t,u,ft,rpar,ipar) 
c                real*8 t,u(n),ft(n),rpar(*)
c                integer ipar(*)
c                ft(1)=...
c
c            rpar,ipar (see below).
c
c   ifdt   = Switch for the computation of df/dt:
c            ifdt = 0: df/dt is computed internally by finite differences,
c                      subroutine fdt is never called (default).
c            ifdt = 1: df/dt is supplied by subroutine fdt.
c
c   solout = Name (external) of user-supplied subroutine exporting the
c            numerical solution during integration. If  iout=1: 
c            solout it is called after every successful step. Supply
c            a dummy subroutine if iout = 0. It must have the form:
c
c               subroutine solout(n,told,tnew,uold,unew,fold,fnew,ucon,
c              1                  intr,rpar,ipar)
c               integer n,intr,ipar(*)
c               real*8 told,tnew,uold(n),unew(n),fold(n),fnew(n),ucon(n)
c               real*8 rpar(*)
c               ...
c               end
c       
c            "intr" serves to interrupt the integration. If solout sets
c             intr .lt. 0, then ROWMAP returns to the calling program.
c
c            solout may produce continuous output by calling the internal
c            subroutine "rowcon" (see below).
c
c   iout   = Gives information on the subroutine solout:
c               iout = 0: subroutine is never called
c               iout = 1: subroutine is used for output (default)
c
c   work   = Array of working space of length "lwork", changed on exit.
c            Serves as working space for all vectors and matrices. 
c            The array is used for optional input. For zero input, the
c            default values are set:
c
c            work(1), work(2) -  Parameters for step size selection. 
c               The new step size is chosen subject to the restriction
c                  work(1) <= hsnew/hsold <= work(2)
c                Default values: work(1)=0.25d0, work(2)=2d0
c   
c            work(3) - The safety factor in step size prediction.
c                Default value: work(3)=0.8d0
c
c            work(4) - UROUND, the machine precision/rounding unit. 
c                Default value: work(4)=1d-16
c
c            work(5) - KTOL, tolerance for the iterative solution of
c                the linear equations. The code keeps the weighted 
c                root-mean-square norm of the residual of the first stage
c                below KTOL/HS, where HS is the current step size
c                (see also ATOL/RTOL).
c                Default value: work(5)=1d-1.
c
c   lwork  = Length of array work in calling program.
c            "lwork" must be at least: 10+n*(mx+11)+mx*(mx+4), 
c            where mx=iwork(3).
c
c   iwork  = Integer working space of length "liwork". The array is used
c            for optional input. For zero input, default values are set:
c            iwork(1)   -  This is the maximal number of allowed steps.
c                Default value is iwork(1)=10000.
c
c            iwork(2)   -  Switch for the choice of integration method:
c                          (see [4], page 110)
c                           1  Method of Shampine
c                           2  Method GRK4T of Kaps-Rentrop 
c                           3  Method of Van Veldhuizen (gamma=1/2)
c                           4  Method of Van Veldhuizen ("D-stable")
c                           5  L-stable Method
c                           6  Method GRK4A of Kaps-Rentrop
c                Default value is iwork(2)=2.
c 
c            iwork(3)   -  The maximum Krylov dimension allowed (=mx). 
c                Default value is iwork(3)=70.
c
c            iwork(4)   -  The logical output unit number "lun" for
c                messages of ROWMAP. Default value is iwork(4)=6.
c                          
c            
c   liwork = Length of array iwork in calling program. 
c            "liwork" must be at least mx+20, where mx is iwork(3). 
c                       
c   rpar, ipar = Real and integer parameters (or parameter arrays)
c            that can be used for communication between your calling
c            program and the subroutines f, fdt, fjac, solout.
c
c   OUTPUT PARAMETERS
c   -----------------
c
c   t     = t-Value where the solution is computed (after successful  
c           return is t=tend)
c
c   u     = Solution at t.
c
c   hs    = Predicted step size from the last accepted step.
c
c   idid  = Reports on success upon return:
c         idid = 1    computation sucessful
c         idid = 2    computation interrupted by solout
c         idid = -1   stepsize too small, i.e. hs < 10*uround*dabs(t)
c         idid = -2   more than iwork(1) steps
c         idid = -3   input is not consistent
c         idid = -4   internal Krylov matrix is repeatedly singular
c
c   iwork(5) = Number of computed (accepted and rejected) steps.
c   iwork(6) = Number of rejected steps.
c   iwork(7) = Number of function evaluations.
c   iwork(8) = Number of jacobian-times-vector products.
c   iwork(9) = Minimum lenght of array work.
c   iwork(10)= Minimum lenght of array iwork.
c
c ROWMAP uses the following basic linear algebra modules (BLAS):
c       Level 1: DAXPY,DCOPY,DNORM2,DROTG,DROT,DDOT,DSCAL
c       Level 2: DGEMV,DTRSV
c
c-----!----------------------------------------------------------------- 
      implicit none
      integer n,lwork,liwork,iwork(liwork)
      integer itol,ijacv,ifdt,iout,ipar(*) ,ifcn,idid
      integer pmq,pmh,pmak1,pmak2,pmak3,pmak4,pmft,pmuu,pmfu0
      integer pmrhs,pml1,pml2,pml3,pml4,pmfm,pmcon,pmscal,pmmax,mx
      real*8 rpar(*) 
      real*8 t,u(n),tend,hs,rtol(*),atol(*),work(lwork),ktol
      external fdt,f,solout,jacv 
c
c Globals.
c
      real*8 uround
      integer cifcn,citol,cijacv,cifdt,ciout,method,lun
      common /rowmap0/uround,cifcn,citol,cijacv,cifdt,ciout,method,lun

c     For statistics.
      integer nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim(4,3)
      common /rowmap1/nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim
c 
c Parameters for step size control.
c
      integer nstpx
      real*8 fac1,fac2,fac3
      common  /rowmap2 /fac1,fac2,fac3,nstpx
c
c Pointers in array "work".
c
      mx=70
      if (iwork(3).gt.19) mx=iwork(3)
      pmq=10
      pmh=pmq+n*mx
      pmak1=pmh+mx*mx
      pmak2=pmak1+n
      pmak3=pmak2+n
      pmak4=pmak3+n
      pmft=pmak4+n
      pmuu=pmft+n
      pmfu0=pmuu+n
      pmrhs=pmfu0+n
      pml1=pmrhs+n
      pml2=pml1+mx
      pml3=pml2+mx
      pml4=pml3+mx
      pmfm=pml4+mx
      pmcon=pmfm+n
      pmscal=pmcon+n
      pmmax=pmscal+n
      iwork(9)=pmmax
c
c Check lwork und liwork.
c
      idid=0
      if (lwork .lt. pmmax) then
          write (lun,9950) lwork, pmmax
 9950     format ('Error in ROWMAP: "lwork"= ',i10,' too small. ', 
     1            'Minimum is',i10,'.')
          idid=-3
      end if
      iwork(10)=mx+20
      if (liwork.lt.iwork(10)) then
          write (lun,9951) liwork,iwork(10)
 9951     format ('Error in ROWMAP: "liwork"=',i6,' too small. ', 
     1           'Minimum is',i6,'.')
          idid=-3
      end if
      if (idid.eq.-3) return
c
c Check optional input.
c
      nstpx=10000
      if (1.lt.iwork(1)) nstpx=iwork(1)
      fac1=work(1)
      fac2=work(2)
      fac3=work(3)
      if (fac1.ge.1.or.fac1.lt.1d-2) fac1=0.25d0
      if (fac2.le.1.or.fac2.gt.1d+2) fac2=2d0
      if (fac3.lt.1d-1.or.fac3.ge.1d0) fac3=0.8d0
c The rounding unit.
      uround=1d-16
      if (work(4).gt.1d-40.and.work(4).lt.1d-5) uround=work(4)
      ktol=1d-1
      if (work(5).gt.0) ktol=work(5)
c     default ROW-method is GRK4T (method=2)
      method=2
      if (1.le.iwork(2).and.6.ge.iwork(2)) method=iwork(2) 
      cifcn=ifcn
      citol=itol
      cijacv=ijacv
      cifdt=ifdt
      ciout=iout
c     logical output unit number for messages
      lun=6
      if (iwork(4).gt.0) lun=iwork(4)
     
c 
c Call to core integrator.
c
      call rowmapc(n,f,t,u,tend,hs,rtol,atol,jacv,fdt,solout,
     1     ktol,work(pmq),work(pmh),
     2     work(pmak1),work(pmak2),
     3     work(pmak3),work(pmak4), work(pmft),
     4     work(pmuu),work(pmfu0),
     5     work(pmrhs),work(pml1),
     6     work(pml2),
     7     work(pml3),work(pml4),
     8     work(pmfm),work(pmcon),work(pmscal),
     9     iwork(20),mx,rpar,ipar,idid) 
c
c Statistics:
c
      iwork(5) = nsteps
      iwork(6) = nstepsr
      iwork(7) = nfeval
      iwork(8) = njacv
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   R O W C O N 
c-----!-----------------------------------------------------------------
c This subroutine is for dense output. It can be called by the
c user-supplied subroutine "solout". A Hermite-interpolation 
c is used to compute an approximation (order 3) "ucon" to the solution
c of the ODE system at time s. The value s should lie in the
c interval [told,tnew].
c 
c Input 
c -----
c  n    = Dimension of the ODE system.
c  s    = Point of evaluation.
c  uold = Numerical solution for the ODE at time told, passed by solout.
c  fold = Value of the right hand side at time told, passed by solout.
c  unew = Numerical solution for the ODE at time tnew, passed by solout.
c  fnew = Value of the right hand side at time tnew, passed by solout.
c  
c Output
c ------
c  ucon = Approximation at time s.
c
      subroutine rowcon(n,s,told,tnew,uold,unew,fold,fnew,ucon)
      implicit none
      integer n
      real*8 s,told,tnew,uold(n),unew(n),fold(n),fnew(n),ucon(n)
      real*8 theta,theta2,thetam1,hs,g
      hs=tnew-told
      theta=(s-told)/hs
      theta2=theta**2d0
      thetam1=theta-1d0
      call dcopy(n,0d0,0,ucon,1)
      g=theta2*(3d0-2d0*theta)
      call daxpy(n,1d0-g,uold,1,ucon,1)
      call daxpy(n,g,unew,1,ucon,1)
      g=hs*theta*thetam1**2d0
      call daxpy(n,g,fold,1,ucon,1)
      g=hs*theta2*thetam1
      call daxpy(n,g,fnew,1,ucon,1)
      return
      end

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E   R O W M A P C
c-----!----------------------------------------------------------------- 
c Here is the core integrator. ROWMAPC prepares the linear stage equations
c and calls subroutine "STAGE" for solving these. ROWMAPC contains the
c step size control. 
c   
      subroutine rowmapc(n,f,t,u,tend,hs,rtol,atol,jacv,fdt,solout,
     1              ktol,q,h,ak1,ak2,ak3,ak4,ft,uu,fu0,rhs,l1,l2,
     2               l3,l4,fm,ucon,scal,kl,mx,rpar,ipar,idid)
      implicit none
      integer mx,n,kl(mx),asv,mk,ipar(*),idid,i,j,intr,info,nsing
      real*8 t,u(n),tend,hs,rtol(*),atol(*),fm(n),ts,ktol,ktol1
      real*8 q(n,mx),h(mx,mx),ak1(n),ak2(n),ak3(n),ak4(n),
     1       ft(n),uu(n), fu0(n),l1(mx),l2(mx),l3(mx),l4(mx),
     2       rhs(n) , dnrm2,unorm,ucon(n),scal(n),
     3       rpar(*),delt,ehg,err,hnew,told
      external f, solout,fdt,jacv

      real*8 a21,a31,a32,c21,c31,c32,c41,c42,c43,b1,b2,b3,b4,e1,e2,e3,
     1   e4,gamma,c2,c3,d1,d2,d3,d4,g21,g31,g32,g41,g42,g43

      logical reached

c
c Common blocks.
c

C     Globals.   
      real*8 uround
      integer ifcn,itol,ijacv,ifdt,iout,method,lun
      common /rowmap0/uround,ifcn,itol,ijacv,ifdt,iout,method,lun
c     For statistics.
      integer nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim(4,3)
      common /rowmap1/nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim

c     For step size control.
      integer nstpx
      real*8 fac1,fac2,fac3
      common  /rowmap2 /fac1,fac2,fac3,nstpx

c     "unorm" is used in subroutine "ROWDQ".
      common /rownorm/unorm
c
c Initializations.
c 
      call rowcoe(method,a21,a31,a32,c21,c31,c32,c41,c42,c43,
     1            b1,b2,b3,b4,e1,e2,e3,e4,gamma,
     2            c2,c3,d1,d2,d3,d4,g21,g31,g32,g41,g42,g43)
      intr=0
      nfeval=0
      nsteps=0
      nstepsr=0
      njacv=0
      nfdt=0
      nsing=0
      reached=.false.
      do 101 i=1,4
         nkrydim(i,1)=mx
         nkrydim(i,2)=0
         nkrydim(i,3)=0
 101  continue
      call f(n,t,u,fm,rpar,ipar)
      nfeval=nfeval+1
      idid=1
c 
c Begin integration, loop for successful steps:
c 
 1    if (t.ge.tend.or.reached) return
      unorm=dnrm2(n,u,1)/dsqrt(dfloat(n))
      if (itol.ne.1) then
c       atol,rtol are scalars
        do 105 i=1,n
           scal(i)=1d0/(atol(1)+rtol(1)*dabs(u(i)))
 105    continue
      else
c       atol,rtol are vectors
        do 106 i=1,n
            scal(i)=1d0/(atol(i)+rtol(i)*dabs(u(i)))
 106    continue
      end if
      if (nsteps.eq.0.and.hs.le.0d0) then
c          Set initial step size.
           hs=dnrm2(n,scal,1)/dsqrt(dfloat(n))
         hs=(1d0/hs)**0.25*1d-1
      end if
      if ((t+hs*1.005).ge.tend) then
        hs=tend-t
        reached=.true.
      end if
c
c Compute the derivative f_t in the nonautonomous case.
c
      if (ifcn.ne.0) then
          nfdt=nfdt+1
          if (ifdt.ne.1) then
c         finite differences:
          delt=dsqrt(uround*max(1d-5,dabs(t)))
          call f(n,t+delt,u,ft,rpar,ipar)
          nfeval=nfeval+1
          call daxpy(n,-1d0,fm,1,ft,1)
          call dscal(n,1d0/delt,ft,1)
      else
c         user-supplied subroutine:
          call fdt(n,t,u,ft,rpar,ipar)
          end if
      end if
c
c Label for rejected steps:
c
 2    ehg=1.D0/(gamma*hs)
      ktol1=ktol/hs
      ts=t
      asv=1
      mk=0
c     Reset the target indices.
      do 102 j=1,mx
        kl(j)=0
  102 continue
      if (nsteps.ge.nstpx) then 
        write (lun,9990) t
        idid=-2
        write (lun,9992) nstpx 
        return
      end if
      nsteps=nsteps+1
c
c 1. stage.
c
      call dcopy (n,fm,1,rhs,1)
      call dscal (n,ehg,rhs,1)
      if (ifcn.ne.0) call daxpy(n,ehg*hs*d1,ft,1,rhs,1)
      call stage(q,h,rhs,l1,kl,asv,mk,mx,n,fm,fu0,u,t,1,ktol1,
     1           ehg,ak1,scal,f,jacv,fdt,rpar,ipar,info)
      if (info.lt.0) goto 3
c
c 2. stage. 
c
      call dcopy(n,u,1,uu,1)
      call daxpy(n,hs*a21,ak1,1,uu,1)
      ts=t+c2*hs 
      call f(n,ts,uu,rhs,rpar,ipar)
      nfeval=nfeval+1
      call daxpy(n,c21,ak1,1,rhs,1)
      call dscal (n,ehg,rhs,1)
      if (ifcn.ne.0) call daxpy(n,ehg*hs*d2,ft,1,rhs,1)
      call stage(q,h,rhs,l2,kl,asv,mk,mx,n,fm,fu0,u,t,2,ktol1,
     1           ehg,ak2,scal,f,jacv,fdt,rpar,ipar,info)
      if (info.lt.0) goto 3
      call daxpy(n,-c21,ak1,1,ak2,1)
c
c 3. stage.
c
      call dcopy(n,u,1,uu,1)
      call daxpy(n,hs*a31,ak1,1,uu,1)
      call daxpy(n,hs*a32,ak2,1,uu,1)
      ts=t+c3*hs 
      call f(n,ts,uu,rhs,rpar,ipar)
      nfeval=nfeval+1
      call dcopy(n,rhs,1,ak4,1)
      call daxpy(n,c31,ak1,1,rhs,1)
      call daxpy(n,c32,ak2,1,rhs,1)
      call dscal(n,ehg,rhs,1)
      if (ifcn.ne.0) call daxpy(n,ehg*hs*d3,ft,1,rhs,1)
      call stage(q,h,rhs,l3,kl,asv,mk,mx,n,fm,fu0,u,t,3,ktol1,
     1           ehg,ak3,scal,f,jacv,fdt,rpar,ipar,info)    
      if (info.lt.0) goto 3
      call daxpy(n,-c31,ak1,1,ak3,1)
      call daxpy(n,-c32,ak2,1,ak3,1)
c
c 4. stage.
c
      call dcopy(n,ak4,1,rhs,1)
      call daxpy(n,c41,ak1,1,rhs,1)
      call daxpy(n,c42,ak2,1,rhs,1)
      call daxpy(n,c43,ak3,1,rhs,1)
      call dscal(n,ehg,rhs,1)
      if (ifcn.ne.0) call daxpy(n,ehg*hs*d4,ft,1,rhs,1)
      call stage(q,h,rhs,l4,kl,asv,mk,mx,n,fm,fu0,u,t,4,ktol1,
     1           ehg,ak4,scal,f,jacv,fdt,rpar,ipar,info)
      if (info.lt.0) goto 3
      nsing=0
      call daxpy(n,-c41,ak1,1,ak4,1)
      call daxpy(n,-c42,ak2,1,ak4,1)
      call daxpy(n,-c43,ak3,1,ak4,1)
c
c New solution: uu = u + sum (h* b.i * ak.i,i=1..4).
c
      call dcopy (n,u,1,uu,1)
      call daxpy (n,hs*b1,ak1,1,uu,1)
      call daxpy (n,hs*b2,ak2,1,uu,1)
      call daxpy (n,hs*b3,ak3,1,uu,1)
      call daxpy (n,hs*b4,ak4,1,uu,1)
c
c Embedded solution: fu0 = sum (hs* e.i * ak.i, i=1..4).
c
      call dcopy (n,0d0,0,fu0,1)
      call daxpy (n,hs*e1,ak1,1,fu0,1)
      call daxpy (n,hs*e2,ak2,1,fu0,1)
      call daxpy (n,hs*e3,ak3,1,fu0,1)
      call daxpy (n,hs*e4,ak4,1,fu0,1)
c
c Error estimate, step size control.
c
       err=0d0
       do 550 i=1,n
          err=err+(fu0(i)*scal(i))**2 
  550  continue
      err=dsqrt(err/n)
      hnew=hs*dmin1(fac2,dmax1(fac1,(1d0/err)**0.25D0 * fac3))
      if (1d-1*hnew.le.dabs(t)*uround) then
          write (lun,9990) t
          idid=-1
          write (lun,9991) 
         return 
      end if
      if (err .lt. 1d0 ) then
c
c Step is accepted.
c
          told=t
          t=t+hs
          call dcopy(n,fm,1,fu0,1)
          call f(n,t,uu,fm,rpar,ipar)
          nfeval=nfeval+1
          if (iout.ne.0) 
     1        call solout(n,told,t,u,uu,fu0,fm,ucon,intr,rpar,ipar)
          call dcopy (n,uu,1,u,1)
          if (intr.lt.0) then
            idid=2
            write (*,9990) t
            write (*,9993)
            return
          end if
          hs=hnew
          goto 1
      else
c
c Step is rejected.
c
          nstepsr=nstepsr+1
          reached=.false.
          hs=hnew
          goto 2
      end if
c
c Matrix is singular.
c
 3    nsing=nsing+1
      write (lun,9995) 
      if (nsing.ge.3) then
        write (lun,9990) ts
        write (lun,9994)
        idid=-5
        return
      end if
      nstepsr=nstepsr+1
      hs=hs/2d0
      goto 2

 9990 format ('Exit of ROWMAP at t=',d10.3)
 9991 format ('Step size too small.')
 9992 format ('More than ',i7,' steps needed.')
 9993 format ('Computation interrupted by "solout".')
 9994 format ('Matrix is repeatedly singular')
 9995 format ('Warning: Matrix is singular.') 
      end           

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   S T A G E 
c-----!-----------------------------------------------------------------
c The subroutine "STAGE" solves the linear equations by using the multiple
c Arnoldi process, see [1,2,3]. The first step incorporates the new right
c hand side into the Krylov subspace. Then the Krylov subspace is extended
c until the residual is small enough or the maximal dimensions are reached.
c Finally the solution "ak" is computed as a linear combination of
c orthogonalized Krylov vectors.
c
c
      subroutine stage(q,h,rhs,l,kl,asv,mk,mx,n,fm,fu0,u,t,st,ktol1,
     1                 ehg,ak,scal,f,jacv,fdt,rpar,ipar,info)
      implicit none
      integer mx,n,mk,i,asv,m,ikrdim,info
      integer kl(mx),st,krydim(4,4),ipar(*)
      logical done
      real*8 q(n,mx),h(mx,mx),rhs(n),l(mx),dnrm2,ktol1,nrmv,nrmrhs
      real*8 fm(n),fu0(n),u(n),t,ehg,dfkt
      real*8 ak(n),scal(n),rpar(*)
      external f,jacv,fdt
c
c     Maximal Krylov dimension increments per stage, see [3]
c     (autonomous/nonautonomous). For the methods shamp, grk4t,
c     veldd, velds and lstab (b3.eq.0)
c
      data (krydim(i,1),i=1,4) /-12,3,1,3/
      data (krydim(i,2),i=1,4) /-16,5,1,5/
c
c    ... and for grk4a (b3.ne.0):
c
      data (krydim(i,3),i=1,4) /-15,3,4,3/
      data (krydim(i,4),i=1,4) /-19,5,4,5/

c     Globals.
      real*8 uround
      integer ifcn,itol,ijacv,ifdt,iout,method,lun
      common /rowmap0/uround,ifcn,itol,ijacv,ifdt,iout,method,lun

c     Statistics.
      integer nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim(4,3)
      common /rowmap1/nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim

      ikrdim=2*(method/6)+ifcn+1
      info=0
      nrmrhs=dnrm2(n,rhs,1)
      if (mk.eq.0) then
          call dcopy(mx,0d0,0,l,1)
      else 
          call orthov(n,mk,q,rhs,l,mx)
      end if
      nrmv=dnrm2(n,rhs,1)
      if (nrmv/(1d-12+nrmrhs) .ge. 1d-8.or.st.lt.4) then
c
c Insert the new (orthogonalized) right hand side "rhs" as 
c new column in Q.
c
         if (mk.gt.0) asv=asv+1
c        Shift last vectors by 1.
         do 100 i=asv-1,1,-1
            call dcopy(n,q(1,mk+i),1,q(1,mk+i+1),1)
 100     continue
         do 101 i=1,mk
            if(kl(i).ge.mk+1) kl(i)=kl(i)+1
 101     continue
         l(mk+1)=nrmv
         call dcopy(n,rhs,1,q(1,mk+1),1)
         call dscal(n,1d0/l(mk+1),q(1,mk+1),1)
      end if
      m=0
      done=.false.
c
c Begin Multiple Arnoldi process.
c
 1000 m=m+1
      if (m.gt.mk) then
          kl(m)=m+asv
      end if
      call kryarn(n,t,u,mk,m,kl,q,h,mx,ehg,l,fm,fu0,f,jacv,rpar,ipar)
      if (m.gt.mk-st+asv.and.m.ge.asv.and.dabs(h(m,m)).lt.uround) then
c        matrix is singular
         info=-1
         return
      end if
      if (m.gt.mk-st+asv.and.m.ge.asv) then
        call krdfkt(n,m,kl,q,h,mx,l,fu0,dfkt,st,asv,ehg,scal)
      else 
        dfkt=1d100
      end if
      if (dfkt.lt.ktol1.and.st.gt.1) done=.true.
      if (dfkt.lt.ktol1.and.m.ge.4) done=.true.
      if(st.gt.1.and.(m-mk).ge.krydim(st,ikrdim)) done=.true.
      if(st.eq.1.and.m.ge.krydim(1,ikrdim)+mx) done=.true.
      if (.not. done) goto 1000
c
c End of multiple Arnoldi process for current stage.
c
      mk=m
c     Used Krylov dimensions.
      nkrydim(st,1)=min0(nkrydim(st,1),mk)
      nkrydim(st,2)=max0(nkrydim(st,2),mk)
      nkrydim(st,3)=nkrydim(st,3)+mk
c     Backsubstitution.
      call dtrsv('u','n','n',m,h,mx,l,1)
c     Compute ak = Q*l.
      call dgemv('n',n,mk,1d0,q,n,l,1,0d0,ak,1)
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   O R T H O V 
c-----!-----------------------------------------------------------------
c This subroutine orthogonalizes v with respect to vectors q(.,i),i=1,m,
c by a modified Gram-Schmitt-Algorithm. Dotproducts are stored in "l",
c i.e. l=Q'v.
c    
c
      subroutine orthov(n,m,q,v,l,mx)
      implicit none
      integer n,m,mx,i
      real*8 q(n,*), v(n),l(mx),s,ddot 
      do 100 i=1,m
          s=ddot(n,q(1,i),1,v,1)
          l(i)=s
          call daxpy(n,-s,q(1,i),1,v,1)
 100  continue
      call dcopy(mx-m,0d0,0,l(m+1),1)
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   R E O R T N 
c-----!-----------------------------------------------------------------
c The vector q(.,m) is incorporated in the Krylov space, it is scaled to
c norm one with h(m,j), where kl(j)=m. The candidates for further subspace
c extentions have to be re-orthogonalized to q(.,m).
c
      subroutine reortn(n,m,kl,q,h,mx)
      implicit none
      integer n,m,kl(*),mx,j0,j
      real*8 q(n,*),h(mx,*),s,ddot
      j0=m
      do 100 j=m-1,1,-1
 100  if (kl(j).ge.m) j0=j
      do 150 j=j0,m-1
          if (kl(j).eq.m.and.h(m,j).gt.1d-14)
     1        call dscal(n,1d0/h(m,j),q(1,m),1)
 150  continue
      do 200 j=j0,m-1
      if (kl(j).gt.m) then
          s=ddot(n,q(1,m),1,q(1,kl(j)),1)
          h(m,j)=s
          call daxpy(n,-s,q(1,m),1,q(1,kl(j)),1)
      end if
 200  continue
      return
      end
c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   K R Y A R N 
c-----!-----------------------------------------------------------------
c This subroutine extends the Krylov subspace by one vector and updates
c the QR-decomposition of matrix "h".
c

      subroutine kryarn(n,t,u,mk,m,kl,q,h,mx,ehg,l,fm,fu0,f,
     1          jacv,rpar,ipar)
      implicit none
      integer n,mk,m,kl(*),mx, i,j,j0

      real*8 uround
      integer ifcn,itol,ijacv,ifdt,iout,method,lun
      common /rowmap0/uround,ifcn,itol,ijacv,ifdt,iout,method,lun

      integer nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim(4,3),ipar(*)
      common /rowmap1/nfeval,nsteps,nstepsr,nfdt,njacv,nkrydim
      real*8  u(n),t,q(n,*),h(mx,*),ehg,l(*), c,s,fm(n),fu0(n)
      real*8  rpar(*)
      external f,jacv
      if (m.gt.mk) call reortn(n,m,kl,q,h,mx)
      l(m) = -l(m)
      j0 = m
      do 50 j=m-1,1,-1
  50    if (kl(j).ge.m) j0=j
c     generate new Givens rotations 
      do 100 j=j0,m-1
        if (m .gt. mk) then
         call drotg(h(j,j),h(m,j),c,s)
         call drot(m-j-1,h(j,j+1),mx,h(m,j+1),mx,c,s)
        else
c        reuse old rotations
         call ztocs(h(m,j),c,s)
        endif
c       apply to right hand side "l"
        call drot(1,l(j),1,l(m),1,c,s)
 100  continue

c
c Multiple Arnoldi-Step with QR-Decomposition.
c
      if (m.le.mk) return
c      Krylov-Step: q(.,kl(m)) = (I - QQ') A q(.,m) 
      if (ijacv.ne.1) then
        call rowdq(n,t,u,q(1,m),fu0,fm,q(1,kl(m)),f,rpar,ipar)
        nfeval=nfeval+1
      else
        call jacv(n,t,u,q(1,m),q(1,kl(m)),rpar,ipar)
       end if
      njacv=njacv+1
      call orthov(n,m,q,q(1,kl(m)),h(1,m),mx)
      h(m,m) = h(m,m)-ehg
c     update the QR-decomposition of matrix "h" 
      do 210 i=2,m
       do 200 j=i-1,1,-1
  200   if (kl(j).ge.i) j0=j
       do 210 j=j0,i-1
         call ztocs(h(i,j),c,s)
         call drot(1,h(j,m),1,h(i,m),1,c,s)
  210 continue
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   K R D F K T 
c-----!-----------------------------------------------------------------
c This subroutine computes the weighted root-mean-square norm of the residual
c
c               r = (I-gh A) Q l - rhs.
c
      subroutine krdfkt(n,m,kl,q,h,mx,r,fu0,dn,st,asv,egh,scal)
      implicit none
      integer n,m,mx,kl(mx),i,j,li,ls,mi,st,asv,j0
      real*8 q(n,*),h(mx,*),r(m),fu0(*),dn,s,qn,dnrm2,egh,gh
      real*8 scal(n),wrmsq

      dn = 0D0
      gh=1d0/egh
      j0 = m
      do 10 j=m,1,-1
  10    if (kl(j).gt.m) j0=j
      ls = m-j0+1

      do 110 i=0,ls-1
       mi = m-i
       s = r(mi)
       do 100 j=0,i-1
  100   s = s-h(mi,m-j)*fu0(ls-j)
       s = s/h(mi,mi)
       fu0(ls-i) = s
       s = dabs(s)
       li = kl(mi)
       if (li.gt.m) then
         qn = dnrm2(n,q(1,li),1)
         h(m+1,mi) = qn
         wrmsq=0d0
         do 101 j=1,n
            wrmsq=wrmsq+(q(j,li)*scal(j))**2 
 101     continue
         wrmsq=dsqrt(wrmsq/n)
         dn = dn + s*wrmsq 
       endif
  110 continue
       dn=dn*gh
      return
      end

c-----!-----------------------------------------------------------------
c --- S U B R O U T I N E   Z T O C S 
c-----!-----------------------------------------------------------------
c Reconstructs old Givens rotations computed by DROTG.
c
      subroutine ztocs(z,c,s)
      implicit none
      real*8 z,c,s
      if (dabs(z).lt.1d0) then
          s=z
          c=dsqrt(1d0-s*s)
      else if (dabs(z).gt.1d0) then
          c=1d0/z
          s=dsqrt(1d0-c*c)
      else 
          c=0d0
          s=1d0
      end if
      return
      end

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E   R O W D Q
c-----!----------------------------------------------------------------- 
c This subroutine approximates the Jacobian-vector product
c by finite differences
c
c     v = (f(u+delta y)-f(u))/delta = A y + O(||u||delta^2)     
c

      subroutine rowdq(n,t,u,y,fu0,fm,v,f,rpar,ipar)
      implicit none
      integer n,ipar(*)
      real*8  t,u(n),y(n),fu0(n),fm(n),v(n),delta
      real*8  rpar(*),eddelta,unorm
      external f
      common /rownorm/ unorm
      delta=1.d-7*dmax1(1d-5,unorm)
      eddelta=1d0/delta
      call dcopy(n,u,1,fu0,1)
      call daxpy(n,delta,y,1,fu0,1)
      call f(n,t,fu0,v,rpar,ipar)
      call daxpy(n,-1d0,fm,1,v,1)
      call dscal(n,eddelta,v,1)
      return
      end

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E   R O W C O E
c-----!----------------------------------------------------------------- 
c This subroutine loads the coefficients of the chosen method.
c
      subroutine rowcoe(method,a21,a31,a32,c21,c31,c32,c41,c42,c43,
     1            b1,b2,b3,b4,e1,e2,e3,e4,gamma,
     2            c2,c3,d1,d2,d3,d4,g21,g31,g32,g41,g42,g43)
      implicit none
      integer method
      real*8 a21,a31,a32,c21,c31,c32,c41,c42,c43,b1,b2,b3,b4,e1,e2,e3,
     1   e4,gamma,c2,c3,d1,d2,d3,d4,g21,g31,g32,g41,g42,g43
      goto (1,2,3,4,5,6),method
 10   continue
      c21=g21/gamma
      c31=g31/gamma
      c32=g32/gamma
      c41=g41/gamma
      c42=g42/gamma
      c43=g43/gamma
      c2=a21
      c3=a31+a32
      d1=gamma
      d2=gamma+g21
      d3=gamma+g31+g32
      d4=gamma+g41+g42+g43
      e1=e1-b1
      e2=e2-b2
      e3=e3-b3
      e4=e4-b4
      return

c -- SHAMP
 1    a21= 1.00000000000000  D0
      a31=0.48000000000000   D0 
      a32=0.12000000000000   D0
      g21=-2.00000000000000  D0
      g31=1.32000000000000   D0 
      g32=0.60000000000000   D0
      g41=-0.05600000000000  D0
      g42= -0.22800000000000 D0 
      g43= -0.10000000000000 D0
      b1=0.29629629629630    D0
      b2=0.12500000000000    D0
      b3=0                   D0
      b4=0.57870370370370    D0
      e1=0                   D0
      e2= -0.04166666666667  D0 
      e3= -0.11574074074074  D0
      e4= 1.15740740740741   D0
      GAMMA= 0.50000000000000D0
      goto 10
c --  GRK4T
 2    g21=-0.27062966775244  D0
      g31=0.31125448329409   D0
      g32=0.00852445628482   D0
      g41=0.28281683204353   D0 
      g42=-0.45795948328073  D0
      g43=-0.11120833333333  D0
      b1=0.21748737165273    D0
      b2=0.48622903799012    D0
      b3=0.00000000000000    D0
      b4=0.29628359035715    D0
      a21=0.46200000000000   D0
      a31=-0.08156681683272  D0
      a32=0.96177515016606   D0
      e1=-0.71708850449933   D0
      e2=1.77617912176104    D0
      e3=-0.05909061726171   D0
      e4=0.00000000000000    D0
      GAMMA=0.23100000000000D0
      goto 10
c -- VELDS
 3    g21=-2.00000000000000  D0
      g31=-1.00000000000000  D0
      g32=-0.25000000000000  D0
      g41=-0.37500000000000  D0
      g42=-0.37500000000000  D0
      g43= 0.50000000000000  D0
      b1=0.16666666666667    D0
      b2=0.16666666666667    D0
      b3=0                   D0
      b4=0.66666666666667    D0
      a21=1.00000000000000   D0
      a31=0.37500000000000   D0
      a32=0.12500000000000   D0
      e1=1.16666666666667    D0
      e2=0.50000000000000    D0
      e3=-0.66666666666667   D0
      e4=0                   D0
      GAMMA=0.50000000000000 D0
      goto 10
c -- VELDD
 4    g21=-0.27170214984937   D0
      g31=0.20011014796684    D0
      g32=0.09194078770500    D0
      g41=0.35990464608231    D0
      g42=-0.52236799086101   D0
      g43=-0.10130100942441   D0
      b1=0.20961757675658     D0
      b2=0.48433148684810     D0
      b3=0.0                  D0
      b4=0.30605093639532     D0
      a21=0.45141622964514    D0
      a31=-0.15773202438639   D0
      a32=1.03332491898823    D0
      e1=-0.74638173030838    D0
      e2=1.78642253324799     D0
      e3=-0.04004080293962    D0
      e4=0.0                  D0
      GAMMA=0.22570811482257  D0
      goto 10
c -- LSTAB
 5    g21=-2.34201389131923   D0
      g31=-0.02735980356646   D0
      g32=0.21380314735851    D0
      g41=-0.25909062216449   D0
      g42=-0.19059462272997   D0
      g43=-0.22803686381559   D0
      b1=0.32453574762832     D0
      b2=0.04908429214667     D0
      b3=0.00000000000000     D0
      b4=0.62637996022502     D0
      a21=1.14564000000000    D0
      a31=0.52092209544722    D0
      a32=0.13429476836837    D0
      e1=0.61994881642181     D0
      e2=0.19268272217757     D0
      e3=0.18736846140061     D0
      e4=0                    D0
      GAMMA=0.57282000000000  D0
      goto 10
c -- GRK4A
 6    a21=0.43800000000000   D0
      a31=0.79692045793846   D0
      a32=0.07307954206154   D0
      g21=-0.76767239548409  D0
      g31=-0.85167532374233  D0
      g32=0.52296728918805   D0
      g41=0.28846310954547   D0
      g42=0.08802142733812   D0
      g43=-0.33738984062673  D0
      b1=0.19929327570063    D0
      b2=0.48264523567374    D0
      b3=0.06806148862563    D0
      b4=0.25000000000000    D0
      e1=0.34632583375795    D0
      e2=0.28569317571228    D0
      e3=0.36798099052978    D0
      e4= 0.0                D0
      GAMMA=0.39500000000000 D0
      goto 10
      end
