c----------------------------------------------------------------------- c Sample driver for the code BIM c----------------------------------------------------------------------- program BIMDO implicit none integer MMAX,MLJACM,lwork,liwork,ltime,lstat parameter(MMAX=1000,MLJACM = 3, & lwork=24 + MMAX*(50+5*MLJACM), & liwork=MMAX+37) double precision y(MMAX),work(lwork) integer iwork(lwork),ijac,iout external fevalr, jevalr, solout character problm*8 parameter(ltime=100,lstat=29) double precision t0,tf,t(0:ltime),h0,h,rtol,atol, & y0(MMAX), rpar(1), scd, mescd,yref(MMAX) integer neqn,mljac,mujac,itol,ndisc,i,j,ierr, idid,ipar(1) integer icount(lstat) integer NSTEPS,NACCEPT,NFAILNEWT,NFAILERR,NITER c----------------------------------------------------------------------- c get the problem dependent parameters c----------------------------------------------------------------------- call prob(problm,neqn,ndisc,t,ijac,mljac,mujac) c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(neqn,y0) c----------------------------------------------------------------------- c read the tolerances and initial stepsize c----------------------------------------------------------------------- write(6,*) 'give the absolute tolerance' read(5,*) atol write(6,*) 'give the relative tolerance' read(5,*) rtol write(6,*) 'give the initial stepsize ' read(5,*) h0 do i=1,neqn y(i)=y0(i) end do do i=1,liwork iwork(i) = 0 end do do i=1,lwork work(i) = 0d0 end do iout = 0 idid = 0 do i=1,lstat icount(i)=0 end do c----------------------------------------------------------------------- c call of the subroutine BIM c----------------------------------------------------------------------- do j=0,ndisc h=h0 t0=t(j) tf=t(j+1) call BIM(neqn,fevalr,t0,tf,y,h,rtol,atol, & jevalr,ijac,mljac,mujac, & work,lwork,iwork,liwork, & rpar,ipar,iout,idid) if (idid.ne.0) then write(6,*) 'ERROR: returned idid =', idid goto 20 endif do i=1,lstat icount(i) = icount(i) + iwork(8+i) end do end do c----------------------------------------------------------------------- c print final solution c----------------------------------------------------------------------- write(6,30) atol,rtol,h0 30 format(//,' we solved the problem with',//, + ' absolute tolerance = ',d10.4,',',/, + ' relative tolerance = ',d10.4,',',/, + ' and initial stepsize = ',d10.4,'.'/) c----------------------------------------------------------------------- c print the integration characteristics c----------------------------------------------------------------------- NSTEPS = 0 NACCEPT = 0 NFAILERR = 0 NITER = 0 DO I=1,5 NITER = NITER + icount(I+4) NSTEPS = NSTEPS + icount(I+9) NACCEPT = NACCEPT + icount(I+14) NFAILERR = NFAILERR + icount(I+19) END DO NFAILNEWT = NSTEPS - NACCEPT - NFAILERR write(6,40) NSTEPS,icount(10),icount(11),icount(12),icount(13), & icount(14), & NACCEPT,NFAILNEWT,NFAILERR, & icount(1),icount(2),icount(3),icount(4), & NITER 40 format( //,'Integration characteristics:',//, + ' # Steps = ',i8,/ + ' # Composition = ',i8,i8,i8,i8,i8,/ + ' # Accept = ',i8,/, + ' # Failnwt = ',i8,/, + ' # Failerr = ',i8,/, + ' # F-eval = ',i8,/, + ' # Jac-eval = ',i8,/, + ' # LU-decomp = ',i8,/, + ' # Linear systems = ',i8,/, + ' # Newt. iterat. = ',i8,///) c----------------------------------------------------------------------- c print error with respect to reference solution c----------------------------------------------------------------------- call solut(neqn,yref) call getscd(mescd,scd,neqn,yref,y,problm,atol,rtol) 20 continue end ccccccccccccccccccccccccccccccc AUXILIARY SUBROUTINES subroutine fevalr(neqn,t,y,dy,ierr,rpar,ipar) double precision t,y,dy,rpar(*) integer neqn,ierr,ipar(*) dimension y(neqn),dy(neqn) call feval(neqn,t,y,dy,ierr,rpar,ipar) return end subroutine jevalr(neqn,t,y,jac,ldjac,ierr,rpar,ipar) double precision t,y,jac,rpar(*) integer neqn,ldjac,ierr,ipar(*) dimension y(neqn),jac(ldjac,neqn) call jeval(neqn,t,y,jac,ldjac,ierr,rpar,ipar) return end subroutine solout(m,t,y,f,k,ord,irtrn) implicit none integer m,k,ord,irtrn double precision t(k),y(m,k),f(m,k) integer i,j 20 format (e22.16,a1) write(20,20) t(k) do i=1,m write(20,20) y(i,k),' ' end do write(20,30) ' ' 30 format(A1,/) irtrn=0 return end c------------------------------------------------------------------------ c The following routines are included in the source file report.f c which is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c auxiliary routines for all drivers c c----------------------------------------------------------------------- subroutine getscd(mescd,scd,neqn,yref,y,problm,atol,rtol) integer neqn double precision mescd, scd, yref(neqn), y(neqn) character*(*) problm logical tolvec double precision rtol, atol c c compute the accuracy of the numerical solution c integer i double precision aerr, aerrmx, aerrl2, rerr, rerrmx, rerrl2 double precision mixed_rerr, mixed_rerrmx logical skipi character*1 skipc integer numa, numr character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) numa = 0 numr = 0 aerrl2 = 0d0 rerrl2 = 0d0 aerrmx = 0d0 rerrmx = 0d0 mixed_rerrmx = 0d0 cwe need the tolerances to computes the mixed error write(*,'('//LF//',a)') 'Numerical solution:' write(*,*) write(*,'('//LF//',a,t41,a)') +' ' , + ' scd ', +' solution component' , + '--------------------------- ignore', +' ' , + 'mixed abs rel abs,rel', +'------------------------------------', +'----- ----- ----- ------- ' do 10 i=1,neqn ConfigureProblem: skipping components in scd computation skipi = + problm.eq.'emep' .and. (i.eq.36 .or. i.eq.38) + .or. problm.eq.'andrews'.and. i.gt.7 + .or. problm.eq.'nand' .and. i.ne.5 + .or. problm.eq.'pump' .and. i.eq.9 + .or. problm.eq.'tba' .and. + (i.ne.49+175 .and. i.ne.130+175 .and. i.ne.148+175) + .or. problm.eq.'fekete' .and. i.gt.60 + .or. problm.eq.'plei' .and. i.gt.14 + .or. problm.eq.'beam' .and. i.gt.40 if (skipi) then skipc = 'y' else skipc = ' ' endif aerr = abs(yref(i)-y(i)) if (.not. skipi) then aerrmx = max(aerrmx,aerr) aerrl2 = aerrl2 + aerr**2 numa = numa + 1 endif if (yref(i).ne.0d0) then rerr = abs((yref(i)-y(i))/(yref(i))) if (.not. skipi) then rerrmx = max(rerrmx,rerr) rerrl2 = rerrl2 + rerr**2 numr = numr + 1 endif endif c mixed relative-absolute error: c each code tries to put it below rtol(i) mixed_rerr = abs((yref(i)-y(i))/ + (atol/rtol+dabs(yref(i)))) mixed_rerrmx = max(mixed_rerrmx,mixed_rerr) if (aerr.eq.0d0) then write(*, + '('//LF// ', ''y('',i3,'') = '','// + 'e24.16e3, t74,a)' + ) i,y(i), skipc elseif (yref(i).eq.0d0) then write(*, + '('//LF// ', ''y('',i3,'') = '','// + 'e24.16e3,t36,f10.2,1x,f10.2, t74,a)' + ) i,y(i),-log10(mixed_rerr),-log10(aerr), skipc else write(*, + '('//LF// ', ''y('',i3,'') = '','// + 'e24.16e3,t36,f10.2,1x,f10.2,1x,f10.2,t74,a)' + ) i,y(i),-log10(mixed_rerr),-log10(aerr),-log10(rerr),skipc endif 10 continue aerrl2 = sqrt(aerrl2) rerrl2 = sqrt(rerrl2) write(*,*) write(*,'('//LF//',a,t36,i10,1x,i10,1x,i10)') + 'used components for scd', neqn,numa, numr write(*,'('//LF//',a,t36,f10.2,1x,f10.2,1x,f10.2)') + 'scd of Y (maximum norm)', + -log10(mixed_rerrmx), -log10(aerrmx), -log10(rerrmx) write(*,'('//LF//',a,t47,f10.2,1x,f10.2)') + 'scd of Y (Euclidian norm)', + -log10(aerrl2), -log10(rerrl2) write(*,*) write(*,*) ConfigureProblem: use relative (default) or absolute error ConfigureProblem: for the scd computation mescd = -log10(mixed_rerrmx) write(*,'('//LF//',a,t36,f10.2)') + 'using mixed error yields mescd', mescd if (problm .eq. 'medakzo') then scd = -log10(aerrmx) write(*,'('//LF//',a,t36,f10.2)') + 'using absolute error yields scd', scd else scd = -log10(rerrmx) write(*,'('//LF//',a,t36,f10.2)') + 'using relative error yields scd', scd endif return end subroutine getfmt(LF, FF, PROMPT) character*(*) LF, FF, PROMPT c c standard Fortran carriage control c LF = '''1''' FF = ''' ''' PROMPT = ' ' c c if supported, you might prefer c c LF = '''''' c FF = '''''' c PROMPT = ',$' c return end integer function lnblnk(string) character*(*) string c c return index of last non blank character in string c integer i do 10 i=len(string),1,-1 if (string(i:i) .ne. ' ') then lnblnk = i return endif 10 continue lnblnk = 0 return end