c-----------------------------------------------------------------------
c     driver for MBDFDAE code, ODE case
c-----------------------------------------------------------------------
      program mbdfdo
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      PARAMETER (ND=500,LWORK=(32+2*ND)*ND+1,LIWORK=ND+3)
      COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,
     +    MAXORD,maxstep,jtotfa
      logical linear
      COMMON /type/ linear
      dimension y(nd),yref(nd),work(lwork),iwork(liwork),mbnd(4),
     +   massbnd(4)
      character driv*9, problm*8, solver*6
      parameter (driv = 'mbdfdo', solver='MEBDF')
c-----------------------------------------------------------------------
c     get the problem dependent parameters
c-----------------------------------------------------------------------
      linear = .false.
      call prob(problm,neqn,tbegin,tend,ijac,mljac,mujac)
      mf = 22
      if (mljac.lt.neqn) mf=24
      if (ijac.eq.1) mf = mf-1
      mbnd(1) = mljac
      mbnd(2) = mujac
      mbnd(3) = mujac+mljac+1
      mbnd(4) = 2*mljac+mujac+1
      massbnd(1) = 0
      maxder  = 7
      itol    = 2
      lout    = 6
      index   = 1
      maxstep = 1000000 
      t = tbegin
c-----------------------------------------------------------------------
c     get the initial values
c-----------------------------------------------------------------------
      call init(neqn,y)
c-----------------------------------------------------------------------
c     read the tolerances and initial stepsize
c-----------------------------------------------------------------------
      write(6,*) 'give relative error tolerance'
      read(5,*)   rtol
      write(6,*) 'give absolute error tolerance'
      read(5,*)   atol
      write(6,*) 'give initial  stepsize '
      read(5,*)   h0
      h = h0

      t0   = tbegin
      tf   = tend
      tout = tend
c-----------------------------------------------------------------------
c     call of the subroutine MBDFDAE
c-----------------------------------------------------------------------
220    continue
       call driver(neqn,t0,h,y,tout,tf,mf,index,lout,lwork,
     &  work,liwork,iwork,mbnd,massbnd,maxder,itol,rtol,atol)
        IF (INDEX.EQ.1) THEN
           INDEX=0
           GOTO 220
        END IF
c-----------------------------------------------------------------------
c     print final solution
c-----------------------------------------------------------------------
      write(6,11)rtol,atol,h0,tf
   11 format(/,' we solved the problem with',//,
     +       ' relative error tolerance = ',d10.4,',',/,
     +       ' absolute error tolerance = ',d10.4,/,
     +       '     and initial stepsize = ',d10.4,//,
     +       ' the numerical solution at t = ',d22.14,' :',/)
      do 20 i=1,neqn
      write(6,21)i,y(i)
   20 continue
   21 format(' y(',i3,') = ',d22.14)
c-----------------------------------------------------------------------
c     print error with respect to reference solution
c-----------------------------------------------------------------------
      call solut(neqn,yref)
      write(6,29)
   29 format(/)
      do 30 i=1,neqn
   30 write(6,31)i,abs(y(i)-yref(i))
   31 format(' absolute error in component ',i3,' = ',d22.14,)
c-----------------------------------------------------------------------
c     print statistics
c-----------------------------------------------------------------------
      write(6,41)NSTEP,NFE,NJE,NDEC
   41 format(/,' # steps =          ',i8,/,
     +         ' # f-eval =         ',i8,/,
     +         ' # Jac-eval =       ',i8,/,
     +         ' # LU-decomp =      ',i8,/)
      end
c-----------------------------------------------------------------------
c     since in MEBDF the format of the subroutines containing the
c     derivative function and its jacobian  differs from the format in
c     the testset, we transform them
c-----------------------------------------------------------------------
      subroutine F(neqn,t,y,dy)
      double precision t,y,dy
      integer neqn
      dimension y(neqn),dy(neqn)
      call feval(neqn,t,y,dy)
      return
      end
      subroutine PDERV(t,y,jac,neqn,ldim)
      double precision t,y,jac
      integer neqn,ldim
      dimension y(neqn),jac(ldim,neqn)
      call jeval(neqn,t,y,jac,ldim)
      return
      end
      subroutine mas(n,am,m)
      dimension am(n,n)
      return
      end