      subroutine limex (n,nzc,nzv,lsc,fcn,jac,t,y,tend,rtol,ymax0,
     1                  hmax,h,ijob,nrw,rw,niw,iw)
c*    begin prologue
      implicit real*8(a-h,o-z)
c
c***********************************************************************
c                                                                      *
c    extrapolation integrator for the solution of                      *
c    linearly implicit differential-algebraic systems of the form      *
c                                                                      *
c                                                                      *
c      (*)   b (t,y) * y' (t) = f (t,y)                                *
c                                                                      *
c            b :(n,n)-matrix,  1.le.rank(b).le.n                       *
c                                                                      *
c***********************************************************************
c
c=======================================================================
c                                                                      =
c     revision 3.2                                                     =
c     latest change: july 30, '87                                      =
c                                                                      =
c=======================================================================
c
c references:
c============
c
c /1/ p. deuflhard, u. nowak:
c     extrapolation integrators for quasilinear implicit ode's
c     university of heidelberg, sfb 123, tech. rep. 332 (1985)
c     in:
c     p. deuflhard, b. enquist (eds):
c     large scale scientific computing
c     birkaeuser, prog.sci.comp. 7, (1987)
c
c /2/ p. deuflhard, e. hairer, j. zugck:
c     one step and extrapolation methods for differential-algebraic
c     systems
c     university of heidelberg, sfb 123, tech. rep. 318 (1985)
c
c=======================================================================
c
c  equation (*) represents an ode, if b is nonsingular, and
c  a differential-algebraic equation (dae), if b is singular.
c  limex uses an iterative realization of the semi-implicit euler-
c  method (eulsim) for the discretization of (*).
c
c
c  external subroutines (to be supplied by the user):
c  --------------------------------------------------
c
c    fcn (n,nzv,t,y,rhs,  right-hand-side f(t,y) of the system, inclu-
c         bv,irv,icv)     ding variable part bv of left-hand-side matrix
c                         b(t,y)
c                         (specify bv(1),irv(1),icv(1) if no variable
c                          entries in the left-hand-side matrix b(t,y)
c                          occur)
c    lsc (nzc,bc,irc,icc) constant part bc of left-hand-side matrix
c                         b(t,y) (non-zero entries only)
c                         (pass a dummy-name, if there are no constant
c                          entries in the left-hand-side matrix b(t,y) )
c    jac (n,m,t,y,yp,a)   analytic jacobian of the residual r(t,y)
c                         r(t,y):= f(t,y)-b(t,y)*y'
c                         (pass a dummy-name, if numerical difference
c                          approximation of the jacobian is wanted)
c      t                  actual position
c      y(n)               values of y at t
c      yp(n)              values of y' at t
c      rhs(n)             values of the right-hand-side functions f(t,y)
c      n                  dimension of the system
c      bv(nzv)            entries of variable part bv of left-hand-side
c                         matrix b(t,y) (sparse mode storage scheme, see
c                         note 1 below)
c      irv(nzv)           integer array containing row-indices of bv
c      icv(nzv)           integer array containing column-indices of bv
c      nzv                number of entries of bv
c      bc(nzc)            entries of constant part bc of left-hand-side
c                         matrix b(t,y) (sparse mode storage scheme, see
c                         note 1 below)
c      irc(nzc)           integer array containing row-indices of bc
c      icc(nzc)           integer array containing column-indices of bc
c      nzc                number of entries of bc
c      a(m,n)             entries of the jacobian a at t,y
c      m                  number of rows of a
c                         (m = n, in the dense matrix case, m < n in
c                          the band matrix case, see note 2 below )
c
c  input parameters (* marks transient parameters):
c  ------------------------------------------------
c
c  * t                  starting point of integration
c  * y(n)               initial values y(1),...,y(n)
c    tend               prescribed final point of integration
c    rtol               prescribed relative precision (.gt.0)
c    ymax0              threshold for relative error control
c                       error is relative (componentwise) with respect
c                       to:       max(ymax(i),ymax0)
c                          where: ymax(i):= maximum of abs(y) for
c                                 component no.i computed so far
c    hmax               maximum permitted stepsize
c                       (hmax=tend-t is recommended)
c  * h                  initial stepsize guess
c                       (for h=0 an initial stepsize guess is internally
c                        generated)
c    ijob(13)           integer vector, controlling the execution of
c                       the job, length:=13
c   * ijob(1)           =0: b is known to be or might be singular
c                       =1: b is known to be nonsingular
c     ijob(2)           =0: numerical difference-approximation of the
c                           jacobian of the right-hand-side f(t,y)
c                           internally generated
c                       =1: analytic jacobian supplied by the user
c     ijob(3)           =0: b-h*a is a dense matrix ( a: jacobian
c                                                     h: stepsize )
c                       =1: b-h*a is a band matrix
c     ijob(4)           = lower bandwidth of b-h*a (excluding the diago-
c                         nal ); set ijob(4)=0, if ijob(3)=0
c     ijob(5)           = upper bandwidth of b-h*a (excluding the diago-
c                         nal ); set ijob(5)=0, if ijob(3)=0
c   * ijob(6)           performance statistics
c                       =0: no output
c                       =1: standard output
c                       =2: additionally integration monitor
c                       =3: additionally enhanced information
c                       =4: additionally iteration monitor ( of itera-
c                           tive realization of discretization )
c     ijob(7)           solution output
c                       =0: no output
c                       =1: initial values and solution values
c                       =2: additionally solution values at intermediate
c                           points chosen by limex
c    nrw                dimension of real work-space,
c                       to be chosen .ge. :
c
c                       (mbh+mb+jm+9)*n+2*nzv+nzc+2*jm*jm+jm
c
c                         mbh = n                  , if b-h*a is dense,
c                             = 2*ijob(4)+ijob(5)+1, if b-h*a is banded,
c                         mb  = n                  , if b-h*a is dense
c                             =   ijob(4)+ijob(5)+1, if b-h*a is banded,
c                         jm  = 5                  , if ijob(1).eq.0
c                             = 7                  , if ijob(1).eq.1
c
c    rw                 real work-space
c    niw                dimension of integer work-space,
c                       to be chosen .ge. :
c
c                       n+2*(nzv+nzc+jm)+km
c
c                         km  = jm-1
c
c    iw                 integer work-space
c
c ==================================================================
c ===  note 1 : example for sparse mode storage scheme in limex  ===
c ==================================================================
c
c                   consider the following case:
c
c                               |0     0    0  |
c                    b(t,y)=    |1     1    0  |
c                               |t+1  -t   y(2)|
c
c    with nzv=3 and nzc=2 storage of bv and bc is done by the following
c    lines:
c
c                   bv (1) = t+1.d0
c                   irv(1) = 3
c                   icv(1) = 1
c                   bv (2) = y(2)
c                   irv(2) = 3
c                   icv(2) = 3
c                   bv (3) = -t
c                   irv(3) = 3
c                   icv(3) = 2
c            c
c                   bc (1) = 1.d0
c                   irc(1) = 2
c                   icc(1) = 1
c                   bc (2) = 1.d0
c                   irc(2) = 2
c                   icc(2) = 2
c
c        (the order of the entries in bv and bc is arbitrary)
c
c ===========================================================
c ===  note 2 : storage of user supplied banded jacobian  ===
c ===========================================================
c
c  in the band matrix case, the following lines may build up the
c  analytic jacobian a;
c  here afl denotes the quadratic matrix a in dense form, and abd the
c  rectangular matrix a in banded form :
c
c                   ml = ijob(4)
c                   mu = ijob(5)
c                   mh = mu+1
c                   do 20 j = 1,n
c                   i1 = max0(1,j-mu)
c                   i2 = min0(n,j+ml)
c                   do 10 i = i1,i2
c                   k = i-j+mh
c                   abd(k,j) = afl(i,j)
c           10      continue
c           20      continue
c
c           the total number of rows needed in  abd  is  ml+mu+1 .
c           the  mu by mu  upper left triangle and the
c           ml by ml  lower right triangle are not referenced.
c
c ===========================================================
c
c
c  output parameters:
c  ------------------
c
c    t                  actual final point of integration
c    y(n)               final values at t
c    h                  stepsize proposal for next integration step
c                       (h.eq.0. ,if limex fails to proceed)
c    ijob(6)    .ge. 0: successful integration
c                       (ijob(6) not altered internally)
c               .lt. 0: error code after a fail run of limex
c                  =-1: more than jrmax stepsize reductions
c                       occurred per basic integration step
c                  =-2: more than nstmax basic integration steps per
c                       interval have been performed
c                  =-3: stepsize proposal for next basic integration
c                       step was too small
c                  =-4: matrix pencil b-h*a is singular:
c                       no or infinitely many solutions exist
c                  =-5: iterative realization of discretization failed
c                       to succeed
c                  =-6: nilpotency of the system is greater than
c                       one
c                  =-7: initial values are inconsistent or nilpotency
c                       of the system is greater than one
c                  =-8: real or integer work-space is exhausted
c                  =-9: the given problem is an algebraic equation
c                       (limex is not suitable in this case)
c
c  output after a successfully completed task only:
c  ------------------------------------------------
c
c   ijob(13)  =   number of performed integration steps
c   ijob( 9)  =   number of function-evaluations (for integration)
c   ijob(12)  =   number of function-evaluations (for jacobian approx.)
c   ijob(10)  =   number of gaussian decompositions
c   ijob(11)  =   number of forward-backward-substitutions
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  the numerical solution of the arising linear equations is done by
c  means of the subroutines lingl and subst ( gauss-algorithm with co-
c  lumn-pivoting and row-interchange ) in the dense matrix case, or by
c  the linpack routines dgbfa and dgbsl in the band matrix case.
c  for special purposes these routines may be substituted.
c
c-----------------------------------------------------------------------
c         this is a driver routine for the core-integrator limex1
c-----------------------------------------------------------------------
c
      double precision y(n),rw(nrw)
      integer iw(niw),ijob(13)
c
      external lsc,fcn,jac
c
      common/limx1/lout
      common/limx2/epmach,small
      common/limx3/nstmax,jrmax,ismax
c
c*    end prologue
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c              initial preparations
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  output unit for integration monitor
c  -----------------------------------
c
      lout=6
c
c  relative machine precision
c  --------------------------
c
      epmach=1.d-16
c
c  square-root of smallest positive machine number
c  -----------------------------------------------
c
      small=1.d-35
c
c  prescribed maximum row number of extrapolation-tableau
c  -----------------------------------------------------
c
      jm=7
      if(ijob(1).eq.1) goto 10
      jm=5
c
c  prescribed maximum column number of extrapolation-tableau
c  --------------------------------------------------------
10    km=jm-1
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  internal parameters
c  (standard values fixed below;
c   to be altered, if necessary, by the skillfull user)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  maximum permitted number of basic integration steps per interval
c  ----------------------------------------------------------------
c
      nstmax=900
c
c  maximum permitted number of stepsize reductions per basic step
c  --------------------------------------------------------------
c  (due to extrapolation-tableau)
c
      jrmax=20
c
c  maximum permitted number of stepsize reductions per basic step
c  --------------------------------------------------------------
c  (due to zero pivot in gaussian decomposition)
c
      ismax=5
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  check for sufficient  work-space and type of left-hand-side
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      ijob(8)=0
      if(nzv.gt.0.and.nzc.gt.0) goto 30
      if(nzv+nzc.eq.0) goto 60
      if(nzc.ne.0) goto 20
      nzc=1
      ijob(8)=2
      goto 30
20    nzv=1
      ijob(8)=1
30    ml=ijob(4)
      mu=ijob(5)
      mb=n
      mbh=n
      if(ijob(3).eq.0) goto 40
      mbh=2*ml+mu+1
      mb=ml+mu+1
40    minrw=(mbh+mb+jm+9)*n+2*nzv+nzc+2*jm*jm+jm+jm*n+n
      miniw=n+2*(nzv+nzc)+2*jm+km
      if(ijob(6).gt.0) write(lout,70) minrw,miniw
      if(minrw.gt.nrw .or. miniw.gt.niw) goto 50
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  build-up of work-space
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      nn=mb*n
      nnh=mbh*n
      njm=n*jm
      jmjm=jm*jm
c
      i2=n+1
      i3=i2+jm
      i4=i3+jm
      i5=i4+km
      i6=i5+nzc
      i7=i6+nzc
      i8=i7+nzv
c
      n2=nnh+1
      n3=n2+nn
      n4=n3+nzc
      n5=n4+nzv
      n6=n5+nzv
      n7=n6+njm
      n8=n7+jmjm
      n9=n8+jmjm
      n10=n9+n
      n11=n10+n
      n12=n11+n
      n13=n12+n
      n14=n13+n
      n15=n14+n
      n16=n15+n
      n17=n16+n
      n18=n17+jm
      n19=n18+njm
c
      iprint=1
      if(ijob(6).eq.0) iprint=0
c
c     call of core integrator limex1
c     ----------------------------------
c
      call limex1(n,mb,mbh,nzc,nzv,jm,km,lsc,fcn,jac,t,y,tend,rtol,ymax0
     1           ,hmax,h,ijob,iw(1),iw(i2),iw(i3),iw(i4),iw(i5),iw(i6),
     2            iw(i7),iw(i8),rw(1),rw(n2),rw(n3),rw(n4),rw(n5),
     3            rw(n6),rw(n7),rw(n8),rw(n9),rw(n10),rw(n11),rw(n12),
     4            rw(n13),rw(n14),rw(n15),rw(n16),rw(n17),rw(n18),
     5            rw(n19))
c
      if(ijob(8).eq.2) nzc=0
      if(ijob(8).eq.1) nzv=0
      if(ijob(6).gt.0) write(lout,100) ijob(13),ijob(9),ijob(12),
     1                               ijob(10),ijob(11)
      if(ijob(6).lt.0.and.iprint.eq.0) write(lout,110) ijob(6)
      return
c
c  fail exit  (work-space exhausted)
c  ---------------------------------
50    continue
      icode=-8
      if(ijob(6).eq.0) write(lout,110) icode
      if(ijob(6).gt.0.and.minrw.gt.nrw) write(lout,80)
      if(ijob(6).gt.0.and.miniw.gt.niw) write(lout,90)
      ijob(6)=icode
      return
60    icode=-9
      if(ijob(6).eq.0) write(lout,110) icode
      if(ijob(6).gt.0) write(lout,120)
      ijob(6)=icode
      return
c
c ----------------------------------------------------------------------
c                  format-statements
c ----------------------------------------------------------------------
c
70    format(/1x,' ********************************',/,
     1        1x,' * minimal required work-space: *',/,
     2        1x,' *   real array     rw(',i6,')  *',/,
     3        1x,' *   integer array  iw(',i6,')  *',/,
     4        1x,' ********************************',/ )
80    format(1x,' *** error *** :  real work-space exhausted',/)
90    format(1x,' *** error *** :  integer work-space exhausted',/)
100   format(/1x,' **********  statistics  **********',/,
     1        1x,' ****    steps     :',i6,'     ****',/,
     2        1x,' ****    f.-ev     :',i6,'     ****',/,
     3        1x,' ****    f.-ev(j)  :',i6,'     ****',/,
     4        1x,' ****    decomp.   :',i6,'     ****',/,
     5        1x,' ****    subst.    :',i6,'     ****',/,
     6        1x,' **********************************',/)
110   format(/1x,' ### error code from limex :',i3,' ###',/)
120   format(/1x,'the given problem is an algebraic equation, limex is n
     1ot suitable',/)
c
c ----------------------------------------------------------------------
c  end of driver routine limex
c ----------------------------------------------------------------------
c
      end
      subroutine limex1(n,mb,mbh,nzc,nzv,jm,km,lsc,fcn,jac,t,y,tend,
     1                  rtol,ymax0,hmax,h,ijob,ipivot,nj,incr,nred,irc,
     2                  icc,irv,icv,b,a,bc,bv0,bvk,dt,d,al,ym,del,
     3                  dz,sm,eta,w1,w2,w3,aj,dtp,yp)
c
      implicit real*8(a-h,o-z)
c
      dimension y(n),b(mbh,n),a(mb,n),bc(nzc),bv0(nzv),bvk(nzv),
     1          dt(n,jm),d(jm,jm),al(jm,jm),ym(n),
     2          del(n),dz(n),sm(n),eta(n),w1(n),w2(n),w3(n),aj(jm),
     3          dtp(n,jm),yp(n)
c
      integer ipivot(n),nj(jm),incr(jm),nred(km),ijob(13)
      integer irc(nzc),icc(nzc),irv(nzv),icv(nzv)
c
      logical lfull,lband,ldir
c
      external lsc,fcn,jac
c
      common/limx4/nsol
      common/limx1/lout
      common/limx2/epmach,small
      common/limx3/nstmax,jrmax,ismax
c
      data epmin/1.d-10/,etadif/1.d-6/,ex1/0.6d0/,ex2/1.5d0/,fmin/1.d-2/
     1     ,four/4.d0/,one/1.d0/,one1/1.01d0/,quart/0.25d0/,rmin/0.9d0/,
     2     ro/0.25d0/,safe/0.5d0/,safein/2.d-2/,ten/1.d1/,
     3     thresh/1.d-1/,zero/0.d0/,half/.5d0/,two/2.d0/,safedm/.8d0/
c
c ----------------------------------------------------------------------
c                        initial preparations
c ----------------------------------------------------------------------
c
      ldir=.false.
      lfull=ijob(3).eq.0
      lband=ijob(3).eq.1
      ml=ijob(4)
      mu=ijob(5)
      ml1=ml+1
c
c
      job=0
      epmach=epmach*ten
      epdiff=dsqrt(epmach)
      hmax=dabs(hmax)
      hmaxu=hmax
      nstep=0
      nfcn=0
      nfcnj=0
      ndec=0
      nsol=0
      h1=tend-t
      init=0
      if(ijob(8).eq.2) goto 1000
      call lsc (nzc,bc,irc,icc)
c
c ----------------------------------------------------------------------
c                          initial scaling
c  (for real life applications to be altered by the skillful user)
c ----------------------------------------------------------------------
c
1000  do 1010 i=1,n
      u=dabs(y(i))
      if(u.lt.ymax0) u=ymax0
      sm(i)=u
1010  eta(i)=etadif
c
c ----------------------------------------------------------------------
c           set parameters for extrapolation and order-control
c ----------------------------------------------------------------------
c
c  stepsize sequence for eulsim
c  ----------------------------
c
      do 1020 j=1,jm
1020  nj(j)=j
c
      do 1030 i=1,n
      yp(i)=epmach
      do 1030 j=1,jm
      dt(i,j)=0.d0
      dtp(i,j)=0.d0
1030  continue
c
      fn=dfloat(n)
      fmb=dfloat(mb)
      fj1=dfloat(nj(1))
      eph=ro*rtol
      aj(1)=fj1
      do 1060 j=2,jm
      j1=j-1
      incr(j1)=0
      nred(j1)=0
      fj=dfloat(nj(j))
      v=aj(j1)+fj-one
      aj(j)=v
      do 1040 k=1,j1
      w=fj/dfloat(nj(k))
1040  d(j,k)=w
      if(j.eq.2) goto 1060
      w=v-fj1
      do 1050 k1=2,j1
      k=k1-1
      u=(aj(k1)-v)/(w*dfloat(k1))
      u=eph**u
1050  al(j1,k)=u
1060  continue
c
c  evaluation of cost coefficients
c  -------------------------------
      costf=one
      costj=fn*costf
      if(lband) costj=fmb*costf
      costs=zero
      costlr=zero
c
      costq=one/quart
      costq=3
      if(costlr.gt.zero .and. costs.gt.zero) costq=costlr/costs
c
      if((costs+costlr+costj).eq.zero) goto 1080
      aj(1)=costj+costlr+(costf+costs)*fj1
      do 1070 j=2,jm
      j1=j-1
1070  aj(j)=aj(j1)+(costf+costs)*(dfloat(nj(j))-one)+costs+costlr
1080  koh=1
      joh=2
1090  continue
      if(joh.ge.jm) goto 1100
      if(aj(joh+1)*one1.gt.aj(joh)*al(joh,koh)) goto 1100
      koh=joh
      joh=joh+1
      goto 1090
1100  k=0
      km=koh
      jmh=joh
      incr(jmh)=-1
      omjo=zero
      if(ijob(6).gt.0) write(lout,9010) rtol,ymax0,km
      if(ijob(7).ge.1) write(lout,9160) (y(i),i=1,n)
      if(ijob(6).ge.3) write(lout,9020)
      if(ijob(6).eq.2) write(lout,9030)
c
c ----------------------------------------------------------------------
c                        basic integration step
c ----------------------------------------------------------------------
c
2000  continue
      if(dabs(h1).le.dabs(t)*epmach*ten) goto 4600
      if(dabs(h1).le.epmach*dabs(h)) go to 4600
      if(ijob(6).eq.2) write(lout,9110) nstep,t
      if(ijob(6).ge.3) write(lout,9100) nstep,nfcn,t,h,k,koh
      if(ijob(7).eq.2.and.nstep.gt.0) write(lout,9000) (y(i),i=1,n)
      if(dabs(h1).ge.one1*dabs(h)) go to 2010
      hr=h
      h=h1
2010  jred=0
      nstc=0
      ising=0
      it=0
      dmh=safedm
      do 2020 k=1,km
2020  incr(k)=incr(k)+1
      hmax=dabs(h1)
      if(hmaxu.lt.hmax) hmax=hmaxu
      call reslim (n,nzv,t,y,dz,ym,bv0,irv,icv,yp,fcn,ijob)
      nfcn=nfcn+1
c
c
c ----------------------------------------------------------------------
c  numerical difference approximation of jacobian  a=-dr/dy (scaled)
c  (feed-back control of discretization and rounding errors)
c ----------------------------------------------------------------------
c
      if(ijob(2).eq.1) goto 2140
      if(lband) goto 2070
c
c   dense jacobian
c   --------------
c
      do 2060 k=1,n
      is=0
2040  w=y(k)
      u=sm(k)*eta(k)
      if(ym(k).gt.zero) u=-u
      y(k)=w+u
      call reslim (n,nzv,t,y,w3,del,bvk,irv,icv,yp,fcn,ijob)
      nfcnj=nfcnj+1
      y(k)=w
      u=sm(k)/u
      sumd=zero
      do 2050 i=1,n
      wz=ym(i)
      wy=del(i)
      w=dabs(wy)
      wy=wz-wy
      wz=dabs(wz)
      if(w.lt.wz) w=wz
      if(w.eq.zero) goto 2050
      w=wy/w
      sumd=sumd+w*w
2050  a(i,k)=wy*u/sm(i)
      sumd=dsqrt(sumd/fn)
      if(sumd.eq.zero .or. is.gt.0) goto 2060
      etad=dsqrt(epdiff/sumd)*eta(k)
      if(etad.gt.fmin) etad=fmin
      if(etad.lt.epmin) etad=epmin
      eta(k)=etad
      is=1
      if(sumd.lt.epmin) goto 2040
2060  continue
      goto 2170
c
c  banded jacobian
c  ---------------
c
2070  do 2080 i=1,mb
      do 2080 k=1,n
2080  a(i,k)=zero
      m=ml+mu+1
      do 2130 jj=1,m
      is=0
2090  nfine=1
      do 2100 k=jj,n,m
      w2(k)=y(k)
      w1(k)=sm(k)*eta(k)
      if(ym(k).gt.zero) w1(k)=-w1(k)
2100  y(k)=w2(k)+w1(k)
      call reslim (n,nzv,t,y,w3,del,bvk,irv,icv,yp,fcn,ijob)
      nfcnj=nfcnj+1
      do 2120 k=jj,n,m
      y(k)=w2(k)
      w1(k)=sm(k)/w1(k)
      sumd=zero
      i1=max0(1,k-mu)
      i2=min0(n,k+ml)
      mh=mu+1-k
      if(mh.lt.0)mh=0
      do 2110 i=i1,i2
      wz=ym(i)
      wy=del(i)
      w=dabs(wy)
      wy=wz-wy
      wz=dabs(wz)
      if(w.lt.wz) w=wz
      if(w.eq.zero) goto 2110
      w=wy/w
      sumd=sumd+w*w
2110  a(mh+i-i1+1,k)=wy*w1(k)/sm(i)
      sumd=dsqrt(sumd/fn)
      if(sumd.eq.zero .or. is.gt.0) goto 2120
      etad=dsqrt(epdiff/sumd)*eta(k)
      if(etad.gt.fmin) etad=fmin
      if(etad.lt.epmin) etad=epmin
      eta(k)=etad
      if(sumd.lt.epmin) nfine=0
2120  continue
      if(nfine.eq.1) goto 2130
      is=1
      goto 2090
2130  continue
      goto 2170
c
c ----------------------------------------------------------------------
c                     analytic expression of jacobian
c ----------------------------------------------------------------------
c
c  subroutine jac should provide jacobian matrix of residual a=dr/dy
c  limex1 uses internally the matrix a=-dr/dy (scaled)
c
2140  call jac (n,mb,t,y,yp,a)
      do 2160 i=1,mb
      do 2150 k=1,n
      kk=i+k-1-mu
      if(lfull) kk=i
      if(kk.lt.1) goto 2150
      if(kk.gt.n) goto 2160
      a(i,k)=-sm(k)*a(i,k)/sm(kk)
2150  continue
2160  continue
c
c
c ----------------------------------------------------------------------
c                     discretization
c ----------------------------------------------------------------------
c
2170  do 2180 k=1,n
2180  dz(k)=dz(k)/sm(k)
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  initial stepsize guess, if h=zero
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      if(init.gt.0.or.h.ne.zero) goto 2500
      anorm1=zero
      do 2200 i=1,n
      row=zero
      do 2190 k=1,n
2190  row=dabs(a(i,k))+row
2200  if(row.gt.anorm1)anorm1=row
      if(anorm1.eq.zero)anorm1=one
      h=safein/anorm1
      if(h.gt.h1)h=h1
      init=1
c
2500  tn=t+h
      ldir=.false.
      if(costq.le.zero) ldir=.true.
      fcm=dabs(h)/hmax
      if(fcm.lt.fmin) fcm=fmin
      dm=dmh*two
      if(ijob(6).ge.4) write(lout,9315)
c
      do 3000 j=1,jmh
      iest=0
      if(j.eq.2) iest=1
      m=nj(j)
      g=h/dfloat(m)
      do 3020 i=1,mb
      do 3010 k=1,n
      b(ml+i,k)=g*a(i,k)
3010  continue
3020  continue
      if(ijob(8).eq.1) goto 3040
      do 3030 i=1,nzv
      ir=irv(i)
      irh=ir
      ic=icv(i)
      if(lband) irh=ml+mu+ir-ic+1
      b(irh,ic)=b(irh,ic)+sm(ic)*bv0(i)/sm(ir)
3030  continue
3040  if(ijob(8).eq.2) goto 3060
      do 3050 i=1,nzc
      ir=irc(i)
      irh=ir
      ic=icc(i)
      if(lband) irh=ml+mu+ir-ic+1
      b(irh,ic)=b(irh,ic)+sm(ic)*bc(i)/sm(ir)
3050  continue
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  semi-implicit euler starting step
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
3060  do 3070 i=1,n
      ym(i)=y(i)
3070  del(i)=g*dz(i)
c
      ifail=0
      if(lfull) call lingl(n,b,ipivot,del,ifail)
      if(lband) call dgbfa(b,mbh,n,ml,mu,ipivot,ifail)
      if(lband) call dgbsl(b,mbh,n,ml,mu,ipivot,del,job)
      ndec=ndec+1
      nsol=nsol+1
      if(ifail.ne.0) goto 4100
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c       check for systems with nilpotency greater than one
c       and/or inconsistent initial values
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      if(ijob(1).eq.1.or.j.gt.2.or.nstep.gt.2) goto 3140
      if(j.gt.1) goto 3110
      delm=zero
      do 3100 i=1,n
      w3(i)=dabs(del(i))
3100  if(w3(i).gt.delm)delm=w3(i)
      goto 3140
3110  do 3120 i=1,n
      if(del(i).eq.zero) goto 3120
      if(dabs(del(i)).lt.small) goto 3120
      w3(i)=w3(i)/dabs(del(i))
      if(dabs(w3(i)).gt.ex2) goto 3120
      if(ijob(6).gt.2.and.nstep.eq.0.and.jred.eq.0) write(lout,9140)i
      if(w3(i)*dabs(del(i)).gt.delm*thresh.and.jred.gt.2) goto 3130
3120  continue
      goto 3140
3130  if(w3(i).lt.ex1) goto 4750
      goto 4760
c
3140  continue
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c  semi-implicit euler discretization (iterative realisation)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      do 3150  i=1,n
3150  ym(i)=ym(i)+del(i)*sm(i)
      if(m.eq.1) goto 3300
      m=m-1
c
c
      do 3299 k=1,m
      th=t+dfloat(k)*g
3155  call fcn (n,nzv,th,ym,del,bvk,irv,icv)
      nfcn=nfcn+1
      if(k.gt.1.or.ldir) goto 3170
      diffb=zero
      if(ijob(8).eq.1) goto 3170
      do 3160 i=1,nzv
      ir=irv(i)
      ic=icv(i)
      diff=bvk(i)-bv0(i)
      diffb=diffb+diff*diff*sm(ic)*sm(ic)/sm(ir)/sm(ir)
3160  continue
      diffb=dsqrt(diffb/dfloat(nzv))
      diffb=dfloat(m)*diffb
3170  do 3180 i=1,n
3180  del(i)=g*del(i)/sm(i)
c
      if(ldir) goto 3200
      if(lfull) call subst (n,b,ipivot,del)
      if(lband) call dgbsl (b,mbh,n,ml,mu,ipivot,del,job)
      nsol=nsol+1
      if(ijob(8).ne.1.and.diffb.gt.rtol*1.d-2) goto 3184
      goto 3290
3184  itfail=0
      do 3188 i=1,n
3188  w1(i)=del(i)
      call iter2 (n,nzv,mbh,koh,b,bv0,bvk,del,sm,w1,
     1             w2,ijob,ipivot,irv,icv,rtol,itfail,
     2             theta,delta,iest,fnh,lfull,lband,
     3             iter,costq)
c
      iest=0
      if(ijob(6).ge.4) write(lout,9321) k,iter,itfail,theta,delta
c
      if(itfail.ge.1) ldir=.true.
      if(itfail.gt.1) goto 3155
      goto 3290
c
3200  continue
      do 3220 i=1,mb
      do 3210 kl=1,n
      b(ml+i,kl)=g*a(i,kl)
3210  continue
3220  continue
      if(ijob(8).eq.1) goto 3240
      do 3230 i=1,nzv
      ir=irv(i)
      irh=ir
      ic=icv(i)
      if(lband) irh=ml+mu+ir-ic+1
      b(irh,ic)=b(irh,ic)+sm(ic)*bvk(i)/sm(ir)
3230  continue
3240  if(ijob(8).eq.2) goto 3260
      do 3250 i=1,nzc
      ir=irc(i)
      irh=ir
      ic=icc(i)
      if(lband) irh=ml+mu+ir-ic+1
      b(irh,ic)=b(irh,ic)+sm(ic)*bc(i)/sm(ir)
3250  continue
3260  continue
c
      ifail=0
      if(lfull) call lingl(n,b,ipivot,del,ifail)
      if(lband) call dgbfa(b,mbh,n,ml,mu,ipivot,ifail)
      if(lband) call dgbsl(b,mbh,n,ml,mu,ipivot,del,job)
      ndec=ndec+1
      nsol=nsol+1
c
      if(ijob(6).ge.4) write(lout,9329) k,ifail
c
      if(ifail.ne.0) goto 4100
c
3290  continue
      iest=0
      do 3295 i=1,n
3295  ym(i)=ym(i)+del(i)*sm(i)
3299  continue
c
c ----------------------------------------------------------------------
c                           extrapolation
c ----------------------------------------------------------------------
c
3300  err=zero
      do 33201 i=1,n
      c=del(i)*sm(i)/g
      v=dtp(i,1)
      dtp(i,1)=c
      if(j.eq.1) goto 33201
      ta=c
      do 33101 k=2,j
      jk=j-k+1
      b1=d(j,jk)
      w=c-v
      u=w/(b1-one)
      c=b1*u
      v=dtp(i,k)
      dtp(i,k)=u
33101 ta=u+ta
      yp(i)=ta
33201 continue
c
3301  continue
      do 3320 i=1,n
      c=ym(i)
      v=dt(i,1)
      dt(i,1)=c
      if(j.eq.1) goto 3320
      ta=c
      do 3310 k=2,j
      jk=j-k+1
      b1=d(j,jk)
      w=c-v
      u=w/(b1-one)
      c=b1*u
      v=dt(i,k)
      dt(i,k)=u
3310  ta=u+ta
      ym(i)=ta
      ta=dabs(ta)
      if(ta.lt.sm(i)) ta=sm(i)
      u=u/ta
      err=err+u*u
3320  continue
      if(j.eq.1) goto 3000
c
c error (scaled root mean square)
c -------------------------------
c
      err=dsqrt(err/fn)
      konv=0
      if(err.le.rtol) konv=1
      err=err/eph
c
c ----------------------------------------------------------------------
c                            order control
c ----------------------------------------------------------------------
c
      k=j-1
      root=one/dfloat(j)
      fc=err**root
      if(fc.lt.fcm) fc=fcm
c
c  optimal order determination
c  ---------------------------
c
      omj=fc*aj(j)
      if(j.gt.2.and.omj*one1.gt.omjo.or.k.gt.joh) goto 3340
      ko=k
      jo=j
      omjo=omj
      fco=fc
3340  continue
      if(j.lt.koh.and.nstep.gt.0) goto 3000
      if(konv.eq.0) goto 3360
      if(ko.lt.k.or.incr(j).lt.0) goto 4500
c
c  possible increase of order
c  --------------------------
c
      if(nred(ko).gt.0) nred(ko)=nred(ko)-1
      fc=fco/al(j,k)
      if(fc.lt.fcm) fc=fcm
      j1=j+1
      if(aj(j1)*fc*one1.gt.omjo) goto 4500
      fco=fc
      ko=jo
      jo=jo+1
      goto 4500
c
c ----------------------------------------------------------------------
c                            convergence monitor
c ----------------------------------------------------------------------
c
3360  red=one/fco
      jk=km
      if(joh.lt.km) jk=joh
      if(k.ge.jk) goto 4000
      if(ko.lt.koh) red=al(koh,ko)*red
      if(al(jk,ko).lt.fco) goto 4000
3000  continue
c
c  stepsize reduction (due to extrapolation tableau)
c --------------------------------------------------
c
4000  red=red*safe
      if(red.ge.rmin) red=rmin
      h=h*red
4010  continue
      if(nstep.eq.0) goto 4020
      nred(koh)=nred(koh)+1
      incr(koh)=-2
4020  jred=jred+1
      if(ijob(6).gt.2) write(lout,9080) jred,red
      if(jred.gt.jrmax) goto 4700
      goto 2500
c
c  stepsize reduction (empirical device)
c  -------------------------------------
c
4100  if(ifail.gt.0.and.ijob(1).eq.0) ising=ising+1
      if(ising.gt.ismax) goto 4730
      if(nstep.eq.0.and.jred.gt.2.and.nstc.gt.2) goto 4760
      hmax=g*fj1*quart
      if(ifail.gt.0) hmax=hmax*safe
      red=hmax/dabs(h)
      h=hmax
      if(jred.gt.0.or.it.gt.0)goto 4020
      goto 4010
c
c ----------------------------------------------------------------------
c              preparations for next basic integration step
c ----------------------------------------------------------------------
c
4500  told=t
      t=tn
      h1=tend-t
      do 4510 i=1,n
      ta=ym(i)
      y(i)=ta
c
c  rescaling
c  ---------
c
      ta=dabs(ta)
      if(ta.lt.sm(i)) goto 4510
      sm(i)=ta
4510  continue
      nstep=nstep+1
      if(nstep.gt.nstmax) go to 4710
c
c
c stepsize prediction
c -------------------
c
      h=h/fco
      koh=ko
      joh=koh+1
      if(dabs(h).gt.dabs(t)*epmach*ten) go to 2000
      go to 4720
c
c ----------------------------------------------------------------------
c                             solution exit
c ----------------------------------------------------------------------
c
4600  h=hr
      hmax=hmaxu
      ijob(13)=nstep
      ijob(9)=nfcn
      ijob(10)=ndec
      ijob(11)=nsol
      ijob(12)=nfcnj
      if(ijob(6).ge.3) write(lout,9100) nstep,nfcn,t,h,k,ko
      if(ijob(6).eq.2) write(lout,9110) nstep,t
      if(ijob(7).gt.0) write(lout,9150) t,(y(i),i=1,n)
      return
c
c ----------------------------------------------------------------------
c                               fail exit
c ----------------------------------------------------------------------
c
4700  continue
      if(ijob(6).gt.0)write(lout,9120)jrmax
      ijob(6)=-1
      goto 4800
4710  continue
      if(ijob(6).gt.0) write(lout,9090) nstmax
      ijob(6)=-2
      goto 4800
4720  continue
      if(ijob(6).gt.0) write(lout,9060)
      ijob(6)=-3
      goto 4800
4730  continue
      if(ijob(6).gt.0) write(lout,9070)
      ijob(6)=-4
      goto 4800
4750  continue
      if(ijob(6).gt.0)write(lout,9050)
      ijob(6)=-6
      goto 4800
4760  continue
      if(ijob(6).gt.0)write(lout,9130)
      ijob(6)=-7
4800  h=zero
      hmax=hmaxu
      return
c
c ----------------------------------------------------------------------
c                  format-statements
c ----------------------------------------------------------------------
c
9000  format(/,1x,' --- current values --- ',/(1x,4d18.9),/)
9010  format(1x,' ------------------------------------------------------
     1----------------',/,
     2       1x,' -limex- :  rel.prec. ',d10.3,'   threshold ',d10.3,'
     3 max. col. ',i3,/,
     4'  ---------------------------------------------------------------
     5-------',/)
9020  format(/,3x,'step',3x,'  f-calls',10x,'t',14x,'h',10x,
     *'used/aimed column',/)
9030  format(/,5x,'step',8x,'t')
9050  format(1x,' problem not solvable with limex, probable reason:',/1x
     *,' nilpotency of the system greater than one',/)
9060  format(/1x,' stepsize reduction failed to succeed '//)
9070  format(/1x,' the matrix-pencil b-h*a is singular, the system is no
     *t solvable',/)
9080  format(/1x,i3,'.reduction, stepsize reduction factor ',d10.3)
9090  format(1x,' more than nstmax=',i3,' integration steps'//)
9100  format(/1x,i5,2x,i9,3x,d15.5,2x,d15.5,i9,i9)
9110  format(/1x,i7,2x,d15.5)
9120  format(1x,' more than jrmax=',i3,' stepsize reductions due to extr
     *apolation tableau')
9130  format(1x,' problem not solvable with limex, probable reason:',/1x
     *,' initial data inconsistent or nilpotency of the system greater t
     *han one',/)
9140  format(1x,' *** warning *** '/1x,' component',i5,'    does not hav
     *e an asymptotic ',/1x,' expansion in the initialization phase '/)
9150  format(/1x,' *** solution at t   =',d18.9,' ***',/(1x,4d18.9))
9160  format(1x,' *** initial values *** ',/(1x,4d18.9))
9315  format(/,2x,'     k  iter  fail',8x,'theta',7x,'delta0')
9321  format(2x,3i6,2d13.3)
9329  format(2x,i6,6x,i6)
c
c ----------------------------------------------------------------------
c                  end of core-integrator limex1
c ----------------------------------------------------------------------
c
      end
      subroutine iter2 (n,nzv,mbh,ko,b,bv0,bvk,del,sm,w1,w2,
     1                  ijob,ipivot,irv,icv,rtol,itfail,
     2                  theta,delta,iest,fnh,lfull,lband,
     3                  iter,costq)
c
      implicit real*8(a-h,o-z)
c
c ----------------------------------------------------------------------
c   subroutine iterat performs the iterative realization of the
c   discretization for variable left-hand-side matrix b(t,y)
c ----------------------------------------------------------------------
c
      dimension b(mbh,n),bv0(nzv),bvk(nzv),del(n),sm(n),w1(n),
     1          w2(n),ijob(13),ipivot(n),irv(nzv),icv(nzv)
c
      logical lfull,lband
c
      common/limx4/nsol
      common/limx1/lout
c
      data zero/0.d0/
      data thtmx1/.25d0/ , thtmx2/.50d0/ , delmx1/.25d0/ , delmx2/.5d0/
c
c  preparations
c  ------------
c
      job=0
      itopt=ifix(sngl(costq))
      itmax=2*itopt
      if(lband) ml=ijob(4)
      if(lband) mu=ijob(5)
      if(lband) ml1=ml+1
      iter=1
      theta=zero
      delta=zero
      if(iest.eq.1) delmax=delmx1
      if(iest.eq.1) thtmax=thtmx1
      if(iest.eq.0) delmax=delmx2
      if(iest.eq.0) thtmax=thtmx2
c
c  descale initial value of iteration
c  ----------------------------------
c
      if(ijob(6).gt.4.and.iest.eq.1) write(lout,150)
      do 10 i=1,n
10    w2(i)=w1(i)*sm(i)
c
c  computation of delta-b*x
c  ------------------------
c
20    do 30 i=1,n
      w1(i)=zero
30    continue
      do 40  i=1,nzv
      ic=icv(i)
      ir=irv(i)
      w1(ir)=w1(ir)+(bv0(i)-bvk(i))*w2(ic)
40    continue
      do 50 i=1,n
50    w1(i)=w1(i)/sm(i)
c
c  solution of the linear equation
c  -------------------------------
c
      if(lfull) call subst(n,b,ipivot,w1)
      if(lband) call dgbsl(b,mbh,n,ml,mu,ipivot,w1,job)
      nsol=nsol+1
c
c  error estimation and calculation of stepsize control parameters
c  ---------------------------------------------------------------
c
      err=zero
      do 80 i=1,n
80    err=err+w1(i)*w1(i)
      err=dsqrt(err/dfloat(n))
      if(iter.gt.2) goto 100
      if(iter.eq.1 .and. err.gt.delmax) goto 220
      if(iter.eq.1) goto 90
      theta=err/delta
      if(theta.gt.thtmax) goto 230
      if(iest.eq.0) goto 100
      zz=theta**dfloat(itopt+2)*delta
      if(zz.gt.rtol*1.d-1) goto 240
      goto 100
90    delta=err
c
100   if(ijob(6).gt.4.and.iest.eq.1) write(lout,160) iter,err
      if(err.le.rtol*1.d-1) goto 130
      iter=iter+1
      do 110 i=1,n
      del(i)=del(i)+w1(i)
110   w2(i)=w1(i)*sm(i)
      if(iter.le.itmax) goto 20
      goto 250
c
c  solution exit
c  -------------
c
130   do 140 i=1,n
140   del(i)=del(i)+w1(i)
      itfail=0
      if(iest.eq.1 .and. iter.gt.itopt) itfail=1
      return
c
c  fail exits
c------------
c
220   itfail=2
      return
c
230   itfail=3
      return
c
240   itfail=4
      return
c
250   itfail=5
      return
c
150   format(10x,' ####       iter           error      ####')
160   format(10x,' ####      ',i4,7x,d12.4,'    ####')
c
c ----------------------------------------------------------------------
c                    end of subroutine iterat
c ----------------------------------------------------------------------
c
      end
      subroutine reslim (neq,nzv,t,u,urhs,ures,bv,ir,ic,up,fcn,ijob)
c
c*********************************************************************
c
c  computes residual
c
c*********************************************************************
c
      implicit real*8(a-h,o-z)
c
c
      dimension u(neq),urhs(neq),ures(neq),up(neq)
      dimension bv(nzv),ir(nzv),ic(nzv)
      dimension ijob(1)
c
      external fcn
c
      call fcn (neq,nzv,t,u,urhs,bv,ir,ic)
c
      do 100 i=1,neq
      ures(i)=urhs(i)
100   continue
c
      if(ijob(8).eq.1) return
c
      do 200 i=1,nzv
      ich=ic(i)
      irh=ir(i)
      ures(irh)=ures(irh)-bv(i)*up(ich)
200   continue
c
      return
      end
      subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
      integer lda,n,ml,mu,ipvt(1),info
      double precision abd(lda,1)
c
c     dgbfa factors a double precision band matrix by elimination.
c
c     dgbfa is usually called by dgbco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c
c     on entry
c
c        abd     double precision(lda, n)
c                contains the matrix in band storage.  the columns
c                of the matrix are stored in the columns of  abd  and
c                the diagonals of the matrix are stored in rows
c                ml+1 through 2*ml+mu+1 of  abd .
c                see the comments below for details.
c
c        lda     integer
c                the leading dimension of the array  abd .
c                lda must be .ge. 2*ml + mu + 1 .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c                0 .le. ml .lt. n .
c
c        mu      integer
c                number of diagonals above the main diagonal.
c                0 .le. mu .lt. n .
c                more efficient if  ml .le. mu .
c     on return
c
c        abd     an upper triangular matrix in band storage and
c                the multipliers which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgbsl will divide by zero if
c                     called.  use  rcond  in dgbco for a reliable
c                     indication of singularity.
c
c     band storage
c
c           if  a  is a band matrix, the following program segment
c           will set up the input.
c
c                   ml = (band width below the diagonal)
c                   mu = (band width above the diagonal)
c                   m = ml + mu + 1
c                   do 20 j = 1, n
c                      i1 = max0(1, j-mu)
c                      i2 = min0(n, j+ml)
c                      do 10 i = i1, i2
c                         k = i - j + m
c                         abd(k,j) = a(i,j)
c                10    continue
c                20 continue
c
c           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
c           in addition, the first  ml  rows in  abd  are used for
c           elements generated during the triangularization.
c           the total number of rows needed in  abd  is  2*ml+mu+1 .
c           the  ml+mu by ml+mu  upper left triangle and the
c           ml by ml  lower right triangle are not referenced.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c     fortran max0,min0
c
c     internal variables
c
      double precision t
      integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
c
c
      m = ml + mu + 1
      info = 0
c
c     zero initial fill-in columns
c
      j0 = mu + 2
      j1 = min0(n,m) - 1
      if (j1 .lt. j0) go to 30
      do 20 jz = j0, j1
         i0 = m + 1 - jz
         do 10 i = i0, ml
            abd(i,jz) = 0.0d0
   10    continue
   20 continue
   30 continue
      jz = j1
      ju = 0
c
c     gaussian elimination with partial pivoting
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 130
      do 120 k = 1, nm1
         kp1 = k + 1
c
c        zero next fill-in column
c
         jz = jz + 1
         if (jz .gt. n) go to 50
         if (ml .lt. 1) go to 50
            do 40 i = 1, ml
               abd(i,jz) = 0.0d0
   40       continue
   50    continue
c
c        find l = pivot index
c
         lm = min0(ml,n-k)
         l = idamax(lm+1,abd(m,k),1) + m - 1
         ipvt(k) = l + k - m
c
c        zero pivot implies this column already triangularized
c
         if (abd(l,k) .eq. 0.0d0) go to 100
c
c           interchange if necessary
c
            if (l .eq. m) go to 60
               t = abd(l,k)
               abd(l,k) = abd(m,k)
               abd(m,k) = t
   60       continue
c
c           compute multipliers
c
            t = -1.0d0/abd(m,k)
            call dscal(lm,t,abd(m+1,k),1)
c
c           row elimination with column indexing
c
            ju = min0(max0(ju,mu+ipvt(k)),n)
            mm = m
            if (ju .lt. kp1) go to 90
            do 80 j = kp1, ju
               l = l - 1
               mm = mm - 1
               t = abd(l,j)
               if (l .eq. mm) go to 70
                  abd(l,j) = abd(mm,j)
                  abd(mm,j) = t
   70          continue
               call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
   80       continue
   90       continue
         go to 110
  100    continue
            info = k
  110    continue
  120 continue
  130 continue
      ipvt(n) = n
      if (abd(m,n) .eq. 0.0d0) info = n
      return
      end
c
      subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
      integer lda,n,ml,mu,ipvt(1),job
      double precision abd(lda,1),b(1)
c
c     dgbsl solves the double precision band system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgbco or dgbfa.
c
c     on entry
c
c        abd     double precision(lda, n)
c                the output from dgbco or dgbfa.
c
c        lda     integer
c                the leading dimension of the array  abd .
c
c        n       integer
c                the order of the original matrix.
c
c        ml      integer
c                number of diagonals below the main diagonal.
c
c        mu      integer
c                number of diagonals above the main diagonal.
c
c        ipvt    integer(n)
c                the pivot vector from dgbco or dgbfa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b , where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgbco has set rcond .gt. 0.0
c        or dgbfa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c     fortran min0
c
c     internal variables
c
      double precision ddot,t
      integer k,kb,l,la,lb,lm,m,nm1
c
      m = mu + ml + 1
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve l*y = b
c
         if (ml .eq. 0) go to 30
         if (nm1 .lt. 1) go to 30
            do 20 k = 1, nm1
               lm = min0(ml,n-k)
               l = ipvt(k)
               t = b(l)
               if (l .eq. k) go to 10
                  b(l) = b(k)
                  b(k) = t
   10          continue
               call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
   20       continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/abd(m,k)
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm,t,abd(la,k),1,b(lb),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm,abd(la,k),1,b(lb),1)
            b(k) = (b(k) - t)/abd(m,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (ml .eq. 0) go to 90
         if (nm1 .lt. 1) go to 90
            do 80 kb = 1, nm1
               k = n - kb
               lm = min0(ml,n-k)
               b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
               l = ipvt(k)
               if (l .eq. k) go to 70
                  t = b(l)
                  b(l) = b(k)
                  b(k) = t
   70          continue
   80       continue
   90    continue
  100 continue
      return
      end
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      subroutine dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision da,dx(1)
      integer i,incx,m,mp1,n,nincx
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
      integer function idamax(n,dx,incx)
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if( n .lt. 1 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      dmax = dabs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
   30 continue
      return
      end
      subroutine lingl(n,a,ipi,x,ifail)
      implicit real*8(a-h,o-z)
c     solution of a system of linear equations, a*x = b, by gaussian
c     elimination with partial pivoting.
      real*8 a(n,n),x(n)
      integer ipi(n)
c
      data  zero/0.d0/
c
c     triangular decomposition
      do 30 k=1,n
      p=zero
      do 40 j=k,n
      if (dabs(a(j,k))-dabs(p)) 40,40,50
 50   p=a(j,k)
      i=j
 40   continue
c
      if(dabs(p))60,60,70
c
c     failure exit
 60   ifail=1
      return
c
 70   r=x(i)
      if (i-k) 80,90,80
 80   x(i)=x(k)
      do 100 j=1,n
      q=a(i,j)
      a(i,j)=a(k,j)
 100  a(k,j)=q
 90   ipi(k)=i
      x(k)=r
      k1=k+1
      if(k1-n) 105,105,30
 105  do 110 i=k1,n
      q=a(i,k)
      if (q.eq.zero) goto 110
      q=q/p
      a(i,k)=q
      x(i)=x(i)-q*r
      if(k1-n) 115,115,110
 115  do 120 j=k1,n
 120  a(i,j)=a(i,j)-q*a(k,j)
 110  continue
  30  continue
c
      k=n
 140  q=x(k)
      if (k.eq.n) goto 125
      k1=k+1
      do 130 i=k1,n
 130  q=q-a(k,i)*x(i)
 125  x(k)=q/a(k,k)
      k=k-1
      if (k) 150,150,140
 150  return
      end
      subroutine subst (n,a,ipi,x)
      implicit real*8(a-h,o-z)
c     solution of a system of linear equations, a*x = b, by gaussian
c     elimination with the matrix in decomposed form. the latter,
c     e.g. is produced by a call of procedure lingl provided the
c     exit ifail is not used
      real*8 a(n,n),x(n)
      integer ipi(n)
      do 1 k=1,n
      i=ipi(k)
      j=k-1
      q=x(i)
      if(i.eq.k) goto 10
      x(i)=x(k)
      x(k)=q
10    if(j.eq.0) goto 1
      do 11 i=1,j
11    q=q-a(k,i)*x(i)
1     x(k)=q
      n1=n+1
      do 2 kk=1,n
      k=n1-kk
      q=x(k)
      if(kk.eq.1) goto 2
      k1=k+1
      do 21 i=k1,n
21    q=q-a(k,i)*x(i)
2     x(k)=q/a(k,k)
      return
      end
C 
C 
C c
C c  example for the call of subroutine limex
C c  integrator for differential-algebraic systems
C c
C       implicit real*8(a-h,o-z)
C       double precision y(3),rw(150)
C       integer iw(30),ijob(13)
C       external fcn,lsc,jac,jdum
C       n=3
C       nzc=1
C       nzv=1
C       t=0.d0
C       y(1)=-1.d0
C       y(2)=5.d0
C       y(3)=4.d0*y(1)/(3.d0*y(1)-8.d0)
C       tend=0.1108d0
C       rtol=1.d-3
C       ymax0=1.d0
C       hmax=tend-t
C       h=1.d-3
C       nrw=150
C       niw=30
C       ijob(1)=0
C       ijob(2)=0
C       ijob(3)=0
C       ijob(4)=0
C       ijob(5)=0
C       ijob(6)=2
C       ijob(7)=1
C       call limex(n,nzc,nzv,lsc,fcn,jac,t,y,tend,rtol,ymax0,hmax,h,ijob,
C      1           nrw,rw,niw,iw)
C       stop
C       end
C c
C       subroutine fcn (n,nzv,t,y,dy,bv,irv,icv)
C       implicit real*8(a-h,o-z)
C       dimension y(n),dy(n),bv(1),irv(1),icv(1)
C       dy(1)=-0.8d0*y(1)+1.d1*y(2)-0.6d0*y(1)*y(3)
C       dy(2)=-1.d1+1.6d0*y(3)/y(2)
C       dy(3)=0.8d0*y(1)+1.6d0*y(3)-0.6d0*y(1)*y(3)
C       bv(1)=1.d0/y(2)
C       irv(1)=2
C       icv(1)=2
C       return
C       end
C c
C       subroutine lsc (nzc,bc,irc,icc)
C       implicit real*8(a-h,o-z)
C       double precision bc(nzc)
C       dimension irc(nzc),icc(nzc)
C       bc(1)=1.d0
C       irc(1)=1
C       icc(1)=1
C       return
C       end
C c
C       subroutine jac (n,mb,t,y,yp,a)
C       implicit real*8(a-h,o-z)
C       double precision y(n),yp(n),a(mb,n)
C       a(1,1)=-0.8d0-0.6d0*y(3)
C       a(1,2)=1.d1
C       a(1,3)=-0.6d0*y(1)
C       a(2,1)=0.d0
C       a(2,2)=(yp(2)-1.6d0*y(3))/y(2)**2
C       a(2,3)=1.6d0/y(2)
C       a(3,1)=0.8d0-0.6d0*y(3)
C       a(3,2)=0.d0
C       a(3,3)=1.6d0-0.6d0*y(1)
C       return
C       end
C       subroutine jdum (n,mb,t,y,yp,a)
C       implicit real*8(a-h,o-z)
C       return
C       end
C 
