c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c generic BiMD 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/bimdd.f c c This is revision c $Id: bimdd.f,v 1.4 2006/10/16 12:42:13 testset Exp $ c c----------------------------------------------------------------------- program BiMDd integer md parameter (md=400) integer lwork, liwork parameter (lwork = 24 + md*(70 + 3*md), + liwork= 40 + md) integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(md), + iwork(liwork),ipar(md+1),idid double precision y(md),dy(md),t(0:100), + h0,rtol(md),atol(md), + work(lwork),rpar(1) logical numjac,nummas,consis,tolvec,rtolvec integer itol,ijac,imas,iout double precision h double precision yref(md) character fullnm*40, problm*8, type*3 character driver*8, solver*8 parameter (driver = 'bimdd', solver='BIMD') external odef,odejac,odemas external daef,daejac,daemas external solout double precision solmax real gettim, timtim, cputim external gettim double precision scd, mescd character fileout*140, filepath*100 character formatout*30, namefile*100 logical printsolout, solref, printout integer nindsol, indsol(md) integer i,j,icount(12:30) do 10 i=1,11 iwork(i) = 0 10 continue do 15 i=1,15 work(i) = 0d0 15 continue iout = 0 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, + neqn,ndisc,t, + numjac,mljac,mujac, + nummas,mlmas,mumas, + ind) if (type.eq.'IDE') then print *, 'BiMDd: ERROR: ', + 'BiMD can not solve IDEs' stop elseif (type.eq.'ODE') then imas = 0 elseif (type.eq.'DAE') then imas = 1 do 20 i=1,neqn if (ind(i).eq.0 .or. ind(i).eq.1) then iwork(9)=iwork(9)+1 elseif (ind(i).eq.2) then iwork(10)=iwork(10)+1 elseif (ind(i).eq.3) then iwork(11)=iwork(11)+1 else print *, 'BiMDd: ERROR: ', + 'BiMD can not solve index ', ind(i), ' problems' stop endif 20 continue do 30 i=2,neqn if (ind(i).lt.ind(i-1)) then print *, 'BiMDd: ERROR: ', + 'BiMD requires the index 1,2,3 variables ', + 'to appear in this order' stop endif 30 continue else print *, 'BiMDd: ERROR: ', + 'unknown Test Set problem type', type stop endif if (numjac) then ijac = 0 else ijac = 1 endif c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(neqn,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(neqn,rtol,atol,tolvec) if (tolvec) then rtolvec = .false. do i = 2,neqn rtolvec = rtolvec .or. + ( rtol(i) .ne. rtol(1)) end do if (rtolvec) then print *, 'BiMDd: ERROR: ', + 'BiMD does not handle vector-valued ', + 'relative error tolerance' stop end if itol = 1 else itol = 0 endif call setoutput(neqn,solref,printsolout, + nindsol,indsol) if (printsolout) then iout = 1 ipar(1) = nindsol do i=1,nindsol ipar(i+1) = indsol(i) end do 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 c----------------------------------------------------------------------- c call of the subroutine BIMD c----------------------------------------------------------------------- do 40 j=12,30 icount(j) = 0 40 continue if (printsolout) then c the initial condition is printed in the oputput file write(formatout,'(a,i5,a)') '(e23.15,',nindsol,'(e23.15))' write(90,formatout) + t(0), (y(indsol(it)),it=1,nindsol) end if timtim = gettim() timtim = gettim() - timtim cputim = gettim() do 60 i=0,ndisc h = h0 if (type.eq.'ODE') then call BIMD(neqn,odef,t(i),t(i+1),y,h,rtol,atol,itol, & odejac,ijac,mljac,mujac, & odemas,imas,mlmas,mumas, & solout,iout, & work,lwork,iwork,liwork, & rpar,ipar,idid) elseif (type.eq.'DAE') then call BIMD(neqn,daef,t(i),t(i+1),y,h,rtol,atol,itol, & daejac,ijac,mljac,mujac, & daemas,imas,mlmas,mumas, & solout,iout, & work,lwork,iwork,liwork, & rpar,ipar,idid) endif if (idid.ne.0) then print *, 'BiMDd: ERROR: ', + 'BiMD returned IDID = ', idid stop endif do 50 j=12,30 icount(j) = icount(j) + iwork(j) 50 continue 60 continue cputim = gettim() - cputim - timtim do 70 j=12,30 iwork(j) = icount(j) 70 continue c----------------------------------------------------------------------- c print numerical solution in endpoint and integration statistics c----------------------------------------------------------------------- printout = .true. if (solref) then call solut(neqn,t(ndisc+1),yref) call getscd(mescd,scd,neqn,yref,y,problm,tolvec,atol,rtol, + printout) else call printsol(neqn,y,problm) end if call report( + driver,problm,solver, + rtol,atol,h0,solmax, + iwork,cputim,scd,mescd +) if (printsolout) then close(90) end if end C======================================================================= C `Test Set for IVP Solvers' solout function for BiMD C======================================================================= subroutine solout(m,k,ord,t0,t,y,f,dd,rpar,ipar,irtrn) implicit none integer m,k,ord,irtrn,ipar(*) double precision t0,t(k),y(m,k),f(m,k),dd(k+1,m),rpar(*) character formatout*30 integer i,nindsol nindsol = ipar(1) write(formatout,'(a,i5,a)') '(e23.15,',nindsol,'(e23.15))' write(90,formatout) + t(k),(y(ipar(i+1),k),i=1,nindsol) irtrn = 0 return end C======================================================================= C `Test Set for IVP Solvers' ODE wrappers for BiMD C======================================================================= c c since in BiMD 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 c----------------------------------------------------------------------- subroutine odef(n,t,y,f,ierr,rpar,ipar) integer n,ierr,ipar(*) double precision t,y(n),f(n),rpar(*) ierr = 0 call feval(n,t,y,y,f,ierr,rpar,ipar) return end subroutine odejac(n,x,y,dfy,ldfy,ierr,rpar,ipar) integer ldfy,n,ierr,ipar(*) double precision x,y(n),dfy(ldfy,n),rpar(*) ierr = 0 call jeval(ldfy,n,x,y,y,dfy,ierr,rpar,ipar) return end subroutine odemas(n,am,lmas,ierr,rpar,ipar) integer n,ierr,lmas,ipar(*) double precision am(lmas,n),rpar(*) return end C======================================================================= C `Test Set for IVP Solvers' DAE wrappers for BiMD C======================================================================= c c since in BiMD the format of the subroutines for the c function f, the Jacobian and the Mass-matrix differ from the c format in the testset, we transform them c c----------------------------------------------------------------------- subroutine daef(n,t,y,f,ierr,rpar,ipar) integer n,ierr,ipar(*) double precision t,y(n),f(n),rpar(*) ierr = 0 call feval(n,t,y,y,f,ierr,rpar,ipar) return end subroutine daejac(n,x,y,dfy,ldfy,ierr,rpar,ipar) integer ldfy,n,ierr,ipar(*) double precision x,y(n),dfy(ldfy,n),rpar(*) ierr = 0 call jeval(ldfy,n,x,y,y,dfy,ierr,rpar,ipar) return end subroutine daemas(n,am,lmas,ierr,rpar,ipar) integer n,lmas,ierr,ipar(*) double precision am(lmas,n),rpar(*) double precision x,y ierr = 0 call meval(lmas,n,x,y,y,am,ierr,rpar,ipar) return end