c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c generic MEBDFDAE driver c c DISCLAIMER: see c http://www.dm.uniba.it/~testset/disclaimer.php c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~testset/src/drivers/mebdfd.f c c This is revision c $Id: mebdfd.f,v 1.6 2006/10/02 10:19:09 testset Exp $ c c----------------------------------------------------------------------- program mebdfd integer maxn parameter (maxn=400) integer lwork, liwork parameter (lwork=(32 + 6*maxn+2 + 2*maxn+1)*maxn+2, + liwork=maxn+14) integer n,ndisc,mbnd(4),masbnd(4),ind(maxn), + iwork(20+4*maxn),mf,maxder,lout,idid double precision y(maxn),dy(maxn),t(0:100),tout,tin, + h0,rtol(maxn),atol(maxn), + work(lwork) logical numjac,nummas,consis,tolvec integer itol,ierr double precision h double precision yref(maxn) character fullnm*40, problm*8, type*3 character driver*8, solver*8 parameter (driver = 'mebdfd', solver='MEBDFDAE') external f, pderv, mas double precision solmax real gettim, timtim, cputim external gettim double precision scd,mescd integer i integer ipar(2) double precision rpar(1) character fileout*140, filepath*100 character formatout*30, namefile*100 logical printsolout, solref, printout integer nindsol, indsol(maxn) idid = 1 ierr = 0 lout = 6 maxder = 7 iwork(14) = 1000000 work(1) = 0d0 c----------------------------------------------------------------------- c check the problem definition interface date c----------------------------------------------------------------------- call chkdat(driver,solver,20060828) c----------------------------------------------------------------------- c get the problem dependent parameters c----------------------------------------------------------------------- call prob(fullnm,problm,type, + n,ndisc,t, + numjac,mbnd(1),mbnd(2), + nummas,masbnd(2),masbnd(3), + ind) if (type.eq.'IDE') then print *, 'MEBDFD: ERROR: ', + 'MEBDFDAE can not solve IDEs' stop elseif (type.eq.'ODE') then masbnd(1) = 0 iwork(1) = n iwork(2) = 0 iwork(3) = 0 elseif (type.eq.'DAE') then masbnd(1) = 1 iwork(1) = 0 iwork(2) = 0 iwork(3) = 0 do 10 i=1,n if (ind(i).eq.0 .or. ind(i).eq.1) then iwork(1) = iwork(1)+1 elseif (ind(i).eq.2) then iwork(2) = iwork(2)+1 elseif (ind(i).eq.3) then iwork(3) = iwork(3)+1 else print *, 'MEBDFD: ERROR: ', + 'MEBDFDAE can not solve index ', ind(i), ' problems' stop endif 10 continue do 20 i=2,n if (ind(i).lt.ind(i-1)) then print *, 'MEBDFD: ERROR: ', + 'MEBDFDAE requires the index 1,2,3 variables ', + 'to appear in this order' stop endif 20 continue else print *, 'MEBDFD: ERROR: ', + 'unknown Test Set problem type', type stop endif if (numjac) then if (mbnd(1).ne.n) then mf = 24 else mf = 22 endif else if (mbnd(1).ne.n) then mf = 23 else mf = 21 endif endif if (mbnd(1).ne.n) then mbnd(3) = mbnd(1) + mbnd(2) + 1 mbnd(4) = 2*mbnd(1) + mbnd(2) + 1 endif if (masbnd(1).eq.0) then masbnd(4) = 0 elseif (masbnd(2).eq.n) then masbnd(4) = n else masbnd(4) = masbnd(2) + masbnd(3) + 1 endif c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(n,t(0),y,dy,consis) c----------------------------------------------------------------------- c read the tolerances and initial stepsize c----------------------------------------------------------------------- call getinp(driver,problm,solver,fullnm, + tolvec,rtol,atol,h0,solmax) call settolerances(n,rtol,atol,tolvec) if (tolvec) then itol = 5 else itol = 2 endif call setoutput(n,solref,printsolout, + nindsol,indsol) if (printsolout) then call getarg(1,filepath) call getarg(2,namefile) if (lnblnk(namefile) .gt. 0) then write(fileout,'(a,a,a,a)') filepath(1:lnblnk(filepath)), + namefile(1:lnblnk(namefile)), solver(1:lnblnk(solver)),'.txt' open(UNIT=90,FILE=fileout) call mfileout(namefile,solver,filepath,nindsol,indsol) else write(fileout,'(a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)), solver(1:lnblnk(solver)),'.txt' open(UNIT=90,FILE=fileout) call mfileout(problm,solver,filepath,nindsol,indsol) end if end if if (printsolout) then write(formatout,'(a,i5,a)') '(e23.15,',nindsol,'(e23.15))' end if c----------------------------------------------------------------------- c call of the subroutine MEBDFDAE c----------------------------------------------------------------------- timtim = gettim() timtim = gettim() - timtim if (printsolout) then c The initial value is printed in the output file write(90,formatout) + t(0), (y(indsol(it)),it=1,nindsol) end if cputim = gettim() do 40 i=0,ndisc h = h0 tin = t(i) 30 continue tout = t(i+1) call mebdf(n,tin,h,y,tout,t(i+1),mf,idid,lout, + lwork,work,liwork,iwork, + mbnd,masbnd,maxder,itol,rtol,atol, + rpar,ipar,f,pderv,mas,ierr) if (ierr.eq.-1) then h = h / 2d0 idid = -1 ierr = 0 goto 30 elseif (idid.eq.1) then if (printsolout) then idid = 3 c ------------------------------------------------- c The approximation at t0+h0 is printed in the output file write(90,formatout) + tin, (y(indsol(it)),it=1,nindsol) c ------------------------------------------------- else idid = 0 end if goto 30 elseif (idid.eq.3) then c ------------------------------------------------- c We print only the point in the considered interval c [t(i) t(i+1)] if ((printsolout).and.(tin.le.tout)) then c ------------------------------------------------- write(90,formatout) + tin, (y(indsol(it)),it=1,nindsol) end if goto 30 elseif (idid.ne.0) then print *, 'MEBDFD: ERROR: ', + 'MEBDFDAE returned IDID = ', idid stop endif c ---------------------------------------------- c The approximation in t(i+1)nis printed in the c output file if (printsolout) then write(90,formatout) + tin, (y(indsol(it)),it=1,nindsol) idid = 3 end if c ---------------------------------------------- 40 continue cputim = gettim() - cputim - timtim c----------------------------------------------------------------------- c print numerical solution in endpoint and integration statistics c----------------------------------------------------------------------- printout = .true. if (solref) then call solut(n,t(ndisc+1),yref) call getscd(mescd,scd,n,yref,y,problm,tolvec,atol,rtol, + printout) else call printsol(n,y,problm) end if call report( + driver,problm,solver, + rtol,atol,h0,solmax, + iwork,cputim,scd,mescd +) end C======================================================================= C `Test Set for IVP Solvers' ODE/DAE wrappers for MEBDFDAE C======================================================================= c c since in MEBDFDAE the format of the subroutines for the c function f and the Jacobian differ from the format c in the testset, we transform them c in the ODE case, subsroutine mas is just a dummy routine c c----------------------------------------------------------------------- subroutine f(n,t,y,ydot,ipar,rpar,ierr) integer n,ipar(*),ierr double precision t,y(n),ydot(n),rpar(*) call feval(n,t,y,y,ydot,ierr,rpar,ipar) return end subroutine pderv(t,y,pd,n,meband,ipar,rpar,ierr) integer meband,n,ipar(*),ierr double precision t,y(n),pd(meband,n),rpar(*) double precision dy call jeval(meband,n,t,y,dy,pd,ierr,rpar,ipar) if (ierr.ne.0) then print *, 'MEBDFD: ERROR: ', + 'can not handle JEVAL IERR' stop endif return end subroutine mas(n,am,masbnd,ipar,rpar,ierr) integer n,masbnd,ipar(*),ierr double precision am(masbnd,n),rpar(*) double precision t,y,dy call meval(masbnd,n,t,y,dy,am,ierr,rpar,ipar) if (ierr.ne.0) then print *, 'MEBDFD: ERROR: ', + 'can not handle MEVAL IERR' stop endif return end