c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c generic MEBDFI 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/mebdfid.f c c This is revision c $Id: mebdfid.f,v 1.7 2006/10/02 10:19:09 testset Exp $ c c----------------------------------------------------------------------- program mebdfid integer maxn parameter (maxn=400) integer lwork, liwork parameter (lwork=(38 + 3*maxn)*maxn+3, + liwork=maxn+14) integer n,ndisc,mbnd(4),masbnd(4),ind(maxn), + iwork(liwork),mf,maxder,lout,idid double precision y(maxn),dy(maxn),t(0:100),tin,tout, + h0,rtol(maxn),atol(maxn), + work(lwork) logical numjac,nummas,consis,tolvec integer ldjac,mljac,mujac integer ldmas,mlmas,mumas double precision mas(maxn*maxn) common /jaccom/ ldjac,mljac,mujac common /mascom/ mas,ldmas,mlmas,mumas integer itol,ierr double precision h double precision yref(maxn) character fullnm*40, problm*8, type*3 character driver*8, solver*8 parameter (driver = 'mebdfid', solver='MEBDFI') external residide, pdervide, residode,pdervode,residdae,pdervdae 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,mlmas,mumas, + ind) mljac=mbnd(1) mujac=mbnd(2) if (type.eq.'IDE') then iwork(1) = 0 iwork(2) = 0 iwork(3) = 0 do 11 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 *, 'MEBDFID: ERROR: ', + 'MEBDFI can not solve index ', ind(i), ' problems' stop endif 11 continue do 21 i=2,n if (ind(i).lt.ind(i-1)) then print *, 'MEBDFID: ERROR: ', + 'MEBDFI requires the index 1,2,3 variables ', + 'to appear in this order' stop endif 21 continue if (mlmas.lt.n) then ldmas = mlmas + mumas + 1 else ldmas = n endif call meval(ldmas,n,t(0),y,dy,mas,ierr,rpar,ipar) c M is constant or dummy for all problems elseif (type.eq.'ODE') then nummas = .false. mlmas = 0 mumas = 0 iwork(1) = n iwork(2) = 0 iwork(3) = 0 elseif (type.eq.'DAE') then nummas = .false. 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 *, 'MEBDFID: ERROR: ', + 'MEBDFI 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 *, 'MEBDFID: ERROR: ', + 'MEBDFI requires the index 1,2,3 variables ', + 'to appear in this order' stop endif 20 continue if (mlmas.lt.n) then ldmas = mlmas + mumas + 1 else ldmas = n endif call meval(ldmas,n,t(0),y,dy,mas,ierr,rpar,ipar) c M is constant or dummy for all problems else print *, 'MEBDFID: 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 else mbnd(4) = n endif ldjac = mbnd(4) 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 MEBDI 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) if (type.eq.'IDE') then call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout, + lwork,work,liwork,iwork, + mbnd,maxder,itol,rtol,atol, + rpar,ipar,pdervide,residide,ierr) elseif (type.eq.'ODE') then call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout, + lwork,work,liwork,iwork, + mbnd,maxder,itol,rtol,atol, + rpar,ipar,pdervode,residode,ierr) elseif (type.eq.'DAE') then call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout, + lwork,work,liwork,iwork, + mbnd,maxder,itol,rtol,atol, + rpar,ipar,pdervdae,residdae,ierr) endif 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 *, 'MEBDFID: ERROR: ', + 'MEBDFI returned IDID = ', idid stop endif c ---------------------------------------------- c The approximation in t(i+1) is 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' IDE wrappers for MEBDFI C======================================================================= c c since in MEBDFI the format of the subroutines for the c function G and its derivatives differ from the format c in the testset, we transform them c c G = f -> residide c c dG/dy +1/con dG/dy' -> pdervide c c----------------------------------------------------------------------- subroutine residide(neqn,t,y,g,dy,ipar,rpar,ierr) integer neqn,ierr,ipar(*) double precision t,y(neqn),dy(neqn),g(neqn),rpar(*) call feval(neqn,t,y,dy,g,ierr,rpar,ipar) return end subroutine pdervide(t,y,pd,n,dy,meband,con,ipar,rpar,ierr) integer meband,n,ipar(*),ierr double precision t,y(n),pd(meband,n),rpar(*),con parameter (maxn=400) double precision dy(n),mas(maxn*maxn) common /jaccom/ ldjac,mljac,mujac common /mascom/ mas,ldmas,mlmas,mumas c compute pd = df/dy ierr = 0 if (mljac.lt.n) then call jeval(ldjac,neqn,t,y,dy,pd(mljac+1,1),ierr,rpar,ipar) else call jeval(ldjac,neqn,t,y,dy,pd(1,1),ierr,rpar,ipar) endif if (ierr.ne.0) then print *, 'MEBDFID: ERROR: ', + 'MEBDFI can not handle JEVAL IERR' stop endif c compute pd = df/dy - 1/con*M if (mljac.lt.n) then c df/dy banded and M banded do 20 j=1,n do 10 i=max(1, mumas+1 + 1-j), + min(mlmas+mumas+1, mumas+1 + neqn-j) pd(mljac+mujac-mumas+i,j) = + pd(mljac+mujac-mumas+i,j) - + 1d0/con*mas((j-1)*ldmas+i) 10 continue 20 continue elseif (mlmas.lt.n) then c df/dy full and M banded do 40 j=1,n do 30 i=max(1, mumas+1 + 1-j), + min(mlmas+mumas+1, mumas+1 + neqn-j) pd(j-1+i-mumas,j) = + pd(j-1+i-mumas,j) - + 1d0/con*mas((j-1)*ldmas+i) 30 continue 40 continue else c df/dy full and M full do 60 j=1,n do 50 i=1,n pd(i,j) = + pd(i,j) - 1d0/con*mas((j-1)*ldmas+i) 50 continue 60 continue endif return end C======================================================================= C `Test Set for IVP Solvers' ODE wrappers for MEBDFI C======================================================================= c c since in MEBDFI the format of the subroutines for the c function G and its derivatives differ from the format c in the testset, we transform them c c G = f - y' -> residode c c dG/dy + c dG/dy' = df/dy - c I -> pdervode c c----------------------------------------------------------------------- subroutine residode(neqn,t,y,g,dy,ipar,rpar,ierr) integer neqn,ierr,ipar(*) double precision t,y(neqn),dy(neqn),g(neqn),rpar(*), d(neqn) call feval(neqn,t,y,d,g,ierr,rpar,ipar) do 10 i=1,neqn g(i) = g(i) - dy(i) 10 continue return end subroutine pdervode(t,y,pd,n,dy,meband,con,ipar,rpar,ierr) integer meband,n,ipar(*),ierr double precision t,y(n),pd(meband,n), rpar(*) double precision dy(n),con integer ldjac,mljac,mujac common /jaccom/ ldjac,mljac,mujac ierr = 0 if (mljac.lt.n) then call jeval(meband,n,t,y,dy,pd(1,1),ierr,rpar,ipar) else call jeval(meband,n,t,y,dy,pd(1,1),ierr,rpar,ipar) endif if (ierr.ne.0) then print *, 'MEBDFID: ERROR: ', + 'MEBDFI can not handle JEVAL IERR' stop endif c compute pd = df/dy - 1d0/con*I if (mljac.lt.n) then do 10 j=1,n pd( mujac+1,j) = + pd( mujac+1,j) - 1d0/con 10 continue else do 20 j=1,n pd(j,j) = + pd(j,j) - 1d0/con 20 continue endif return end C======================================================================= C `Test Set for IVP Solvers' DAE wrappers for MEBDFI C======================================================================= c c since in MEBDFI the format of the subroutines for the c function G and its derivatives differ from the format c in the testset, we transform them c c G = f - My' -> residdae c c dG/dy + c dG/dy' = df/dy - cM -> pdervdae c c----------------------------------------------------------------------- subroutine residdae(neqn,t,y,g,dy,ipar,rpar,ierr) integer neqn,ierr,ipar(*) double precision t,y(neqn),dy(neqn),g(neqn),rpar(*) parameter (maxn=400) double precision mas(maxn*maxn) common /mascom/ mas,ldmas,mlmas,mumas integer i,j c compute f in g call feval(neqn,t,y,dy,g,ierr,rpar,ipar) c compute g := f - M*y' if (mlmas.lt.neqn) then do 20 j=1,neqn do 10 i=max(1,j-mumas),min(neqn,j+mlmas) g(i) = + g(i) - mas((j-1)*ldmas+mumas+1-j+i)*dy(j) 10 continue 20 continue else do 40 j=1,neqn do 30 i=1,neqn g(i) = + g(i) - mas((j-1)*ldmas+i)*dy(j) 30 continue 40 continue endif return end subroutine pdervdae(t,y,pd,n,dy,meband,con,ipar,rpar,ierr) integer meband,n,ipar(*),ierr double precision t,y(n),pd(meband,n),rpar(*), con parameter (maxn=400) double precision dy(n),mas(maxn*maxn) common /jaccom/ ldjac,mljac,mujac common /mascom/ mas,ldmas,mlmas,mumas integer i,j c compute pd = df/dy ierr = 0 if (mljac.lt.n) then call jeval(ldjac,n,t,y,dy,pd(1,1),ierr,rpar,ipar) else call jeval(ldjac,n,t,y,dy,pd(1,1),ierr,rpar,ipar) endif if (ierr.ne.0) then print *, 'MEBDFI: ERROR: ', + 'MEBDFI can not handle JEVAL IERR' stop endif c compute pd = df/dy - 1/con*M if (mljac.lt.n) then c df/dy banded and M banded do 20 j=1,n do 10 i=max(1, mumas+1 + 1-j), + min(mlmas+mumas+1, mumas+1 + n-j) pd(mujac-mumas+i,j) = + pd(mujac-mumas+i,j) - + 1d0/con*mas((j-1)*ldmas+i) 10 continue 20 continue elseif (mlmas.lt.n) then c df/dy full and M banded do 40 j=1,n do 30 i=max(1, mumas+1 + 1-j), + min(mlmas+mumas+1, mumas+1 + n-j) pd(j-1+i-mumas,j) = + pd(j-1+i-mumas,j) - + 1d0/con*mas((j-1)*ldmas+i) 30 continue 40 continue else c df/dy full and M full do 60 j=1,n do 50 i=1,n pd(i,j) = + pd(i,j) - 1d0/con*mas((j-1)*ldmas+i) 50 continue 60 continue endif return end