c c The original file report.f has been modified in order to include c the solver GAM c----------------------------------------------------------------------- c c This file 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 DISCLAIMER: see c http://www.dm.uniba.it/~testset/disclaimer.html c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~testset/src/auxil/report.f c c This is revision c $Id: report.F,v 1.8 2003/09/04 10:12:35 testset Exp $ c c----------------------------------------------------------------------- c c This source file contains the common subprograms used by c all Test Set driver programs. c c To get CPU times you should configure gettim: c use your favorite editor to search for: ConfigureTiming c this is the only thing you need to configure. c c However, if you wish to add problems / solvers / drivers c you have to configure more: c in that case use your favorite editor to search for: c c ConfigureProblem c ConfigureSolver c ConfigureDriver c c----------------------------------------------------------------------- ConfigureTiming: let gettim return CPU time real function gettim() c c gettim is supposed to deliver the cpu time (user time) c beware of the clock resolution: typically 1/50 - 1/100 s c c in Fortran 77 timing is highly machine/compiler specific c on some machines you might even have to use a C function instead c c dummy; returns 0 write(*,*) 'WARNING: report.f: ', + 'currently no timing; activate timing by configuring', + 'the function gettim in the file report.f' gettim = 0 c best bet for UNIX c c real*4 etime c real*4 total, tarray(2) c total = etime(tarray) c gettim = tarray(1) c Cray c c real second c gettim = second() return end subroutine chkdat(driver,solver,drivd) character*(*) driver, solver integer drivd c c check date of problem interface c c to be full proof we should check all release dates c i.e. of the driver, solver, problem, and of this file: report.f c c however, the most important date: the release date of the solver c is not available in the solver c so we only do a sanity check on the problem interface date c integer lnblnk, pidate integer probd probd=pidate() probd=probd if (probd.ne.drivd) then write(*,*) 'ERROR: report.f: ', + 'unknown Test Set problem interface date', probd if (probd.gt.drivd) then write(*,*) 'you should probably get and ', + 'install a more recent ', + driver(:lnblnk(driver)),'.f' else write(*,*) 'you should probably get and ', + 'install a more recent problem' endif write(*,*) 'from http://www.dm.uniba.it/~testset/' stop endif return end subroutine getinp(driver,problm,solver,fullnm, + tolvec,rtol,atol,h0,solmax) character*(*) driver, problm, solver, fullnm double precision rtol(*), atol(*), h0, solmax logical tolvec c c get input c check whether driver, problem, and solver are supported c read rtol, atol, h0, solmax c integer i double precision tol logical error logical batch integer lnblnk external batch,lnblnk character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) error = .false. ConfigureProblem: straightforward if (.not.( + problm.eq.'chemakzo' .or. + problm.eq.'hires' .or. + problm.eq.'vdpol' .or. + problm.eq.'vdpolm' .or. + problm.eq.'rober' .or. + problm.eq.'orego' .or. + problm.eq.'e5' .or. + problm.eq.'pollu' .or. + problm.eq.'ringmod' .or. + problm.eq.'andrews' .or. + problm.eq.'transamp' .or. + problm.eq.'medakzo' .or. + problm.eq.'emep' .or. + problm.eq.'nand' .or. + problm.eq.'pump' .or. + problm.eq.'wheel' .or. + problm.eq.'tba' .or. + problm.eq.'caraxis' .or. + problm.eq.'fekete' .or. + problm.eq.'plei' .or. + problm.eq.'water' .or. + problm.eq.'crank' .or. + problm.eq.'beam' + )) +then write(*,*) 'ERROR: report.f: unknown testset problem: ', problm error = .true. endif ConfigureDriver: straightforward if (.not.( + driver.eq.'dassld' .or. + driver.eq.'dopri5d' .or. + driver.eq.'gamd' .or. + driver.eq.'gamdd' .or. + driver.eq.'mebdfd' .or. + driver.eq.'mebdfid' .or. + driver.eq.'psided' .or. + driver.eq.'psoded' .or. + driver.eq.'radaud' .or. + driver.eq.'radau5d' .or. + driver.eq.'voded' + )) +then write(*,*) 'ERROR: report.f: unknown testset driver: ', driver error = .true. endif ConfigureSolver: straightforward if (.not.( + solver.eq.'DASSL' .or. + solver.eq.'DOPRI5' .or. + solver.eq.'GAM' .or. + solver.eq.'GAMD' .or. + solver.eq.'MEBDFDAE' .or. + solver.eq.'MEBDFI' .or. + solver.eq.'PSIDE' .or. + solver.eq.'PSODE' .or. + solver.eq.'RADAU' .or. + solver.eq.'RADAU5' .or. + solver.eq.'VODE' + )) +then write(*,*) 'ERROR: report.f: unknown testset solver: ', solver error = .true. endif if (error) stop write(*,'('//LF//'a,//,'//LF//'a,1x,a,1x,a,1x,a,/)') + ' Test Set for IVP Solvers (release 2.2)', + ' Solving', + fullnm(:lnblnk(fullnm)), + 'using', + solver(:lnblnk(solver)) tolvec = .false. rtol(1) = 0d0 atol(1) = 0d0 h0 = 0d0 solmax = 0d0 if (.not. batch()) then ConfigureSolver: read rtol/atol ConfigureProblem: read special tolerance write(*,'('//LF//'a)') 'User input:' write(*,*) if (problm .eq. 'pump') then write(*,'('//LF//'a'//PROMPT//')') + 'give (pump) error tolerance: ' read *, tol rtol(1) = tol atol(1) = tol elseif (solver .eq. 'PSODE') then write(*,'('//LF//'a'//PROMPT//')') + 'give (mixed) error tolerance: ' read *, tol rtol(1) = tol atol(1) = tol else write(*,'('//LF//'a'//PROMPT//')') + 'give relative error tolerance: ' read *, rtol(1) write(*,'('//LF//'a'//PROMPT//')') + 'give absolute error tolerance: ' read *, atol(1) endif ConfigureSolver: read h0, only if used by solver if ( + solver.eq.'GAM' .or. + solver.eq.'GAMD' .or. + solver.eq.'MEBDFDAE' .or. + solver.eq.'MEBDFI' .or. + solver.eq.'PSODE' .or. + solver.eq.'RADAU' .or. + solver.eq.'RADAU5' + ) + then write(*,'('//LF//'a'//PROMPT//')') + 'give initial stepsize: ' read *, h0 endif ConfigureSolver: read solmax, only if used by solver if (solver.eq.'PSODE') then write(*,'('//LF//'a'//PROMPT//')') + 'give solmax: ' read *, solmax endif else ConfigureSolver: generally CWI only (i.e. normally batch() .eq. .false.) ConfigureSolver: read rtol, atol, h0, solmax; same as above if (solver .eq. 'PSODE') then read *, rtol(1), atol(1), h0, solmax elseif ( + solver.eq.'GAM' .or. + solver.eq.'GAMD' .or. + solver.eq.'MEBDFDAE' .or. + solver.eq.'MEBDFI' .or. + solver.eq.'RADAU' .or. + solver.eq.'RADAU5' + ) then read *, rtol(1), atol(1), h0 elseif (problm .eq. 'pump') then read *, tol rtol(1) = tol atol(1) = tol else read *, rtol(1), atol(1) endif endif c a hack to allow input of e.g. rtol = 10d0 ** (-1.5d0) exactly c negative input is assumed to be the log10 value if (atol(1).lt.0d0) atol(1)=10d0**atol(1) if (rtol(1).lt.0d0) rtol(1)=10d0**rtol(1) if (h0 .lt.0d0) h0 =10d0**h0 if (solmax .lt.0d0) solmax =10d0**solmax if (batch()) then ConfigureSolver: generally CWI only (i.e. normally batch() .eq. .false.) ConfigureSolver: print rtol, atol, h0, solmax; see read above if (problm .eq. 'pump') then write(*,'('//LF//'a,e12.3)') + 'with (pump) error tolerance: ', rtol(1) elseif (solver .eq. 'PSODE') then write(*,'('//LF//'a,e12.3)') + 'with (mixed) error tolerance: ', rtol(1) else write(*,'('//LF//'a,e12.3)') + 'with relative error tolerance: ', rtol(1) write(*,'('//LF//'a,e12.3)') + ' absolute error tolerance: ', atol(1) endif if ( + solver.eq.'GAM' .or. + solver.eq.'GAMD' .or. + solver.eq.'MEBDFDAE' .or. + solver.eq.'MEBDFI' .or. + solver.eq.'PSODE' .or. + solver.eq.'RADAU' .or. + solver.eq.'RADAU5' + ) then write(*,'('//LF//'a,e12.3)') + ' initial stepsize: ', h0 endif if (solver.eq.'PSODE') then write(*,'('//LF//'a,e12.3)') + ' solmax: ', + solmax endif endif ConfigureProblem: special (vector) tolerances if (problm.eq.'pump') then do 10 i=9,6,-1 atol(i) = atol(1) rtol(i) = rtol(1) 10 continue c c charges are much smaller: c do 20 i=5,1,-1 atol(i) = 1d-6*atol(1) rtol(i) = rtol(1) 20 continue tolvec = .true. elseif (problm.eq.'water') then do 30 i=2,49 atol(i) = atol(1) rtol(i) = rtol(1) 30 continue c c pressures are much larger: c do 40 i=37,49 atol(i) = 1d6*atol(1) 40 continue tolvec = .true. endif write(*,*) return end subroutine getscd(mescd,scd,neqn,yref,y,problm,tolvec,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 ned_rerrmx = 0d0 cwe need the tolerances to computes the mixed error if (.not.tolvec) then do i=2,neqn atol(i)=atol(1) rtol(i)=rtol(1) enddo endif 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(i)/rtol(i))+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) c write(*,'('//LF//'a,t36,f10.2,1x,f10.2)') c + 'scd of Y (Euclidian norm)', c + -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,t47,f10.2)') + 'using absolute error yields scd', scd else scd = -log10(rerrmx) write(*,'('//LF//'a,t58,f10.2)') + 'using relative error yields scd', scd endif return end subroutine report( + driver, problm, solver, + rtol, atol, h0, solmax, + iwork, cputim, scd,mescd +) character*(*) driver, problm, solver integer iwork(*) double precision rtol, atol, h0, solmax real cputim double precision scd,mescd c c print integration characteristics including CPU time c integer nsteps, naccpt, nfe, njac, nlu logical batch external batch character*(*) FS, RS c c for LaTeX tabular c c parameter (FS = '''&''', RS = '''\\\\''') c parameter (FS = ''' ''', RS = '''''') character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) ConfigureSolver: get integration characteristics from iwork if (solver.eq.'DASSL') then nsteps = iwork(11)+iwork(14)+iwork(15) naccpt = iwork(11) nfe = iwork(12) njac = iwork(13) nlu = 0 elseif (solver.eq.'DOPRI5') then nsteps = iwork(18) naccpt = iwork(19) nfe = iwork(17) njac = 0 nlu = 0 elseif (solver.eq.'GAMD' .or. solver.eq.'GAM' ) then nsteps = 0 do 12 i=12,23 nsteps = nsteps + iwork(i) 12 continue naccpt = iwork(12)+iwork(13)+iwork(14)+iwork(15) nfe = iwork(10) njac = iwork(11) nlu = iwork(24) elseif (solver.eq.'MEBDFDAE') then nsteps = iwork(5)+iwork(6) naccpt = iwork(5) nfe = iwork(7) njac = iwork(8) nlu = iwork(9) elseif (solver.eq.'MEBDFI') then nsteps = iwork(5)+iwork(6) naccpt = iwork(5) nfe = iwork(7) njac = iwork(8) nlu = iwork(9) elseif (solver.eq.'PSIDE') then nsteps = iwork(15) naccpt = iwork(15)-iwork(16)-iwork(17)-iwork(18)-iwork(19) nfe = iwork(11) njac = iwork(12) nlu = iwork(13) elseif (solver.eq.'PSODE') then nsteps = iwork(1) naccpt = iwork(2) nfe = iwork(3) njac = iwork(4) nlu = iwork(5) elseif (solver.eq.'RADAU' .or. solver.eq.'RADAU5') then nsteps = iwork(16) naccpt = iwork(17) nfe = iwork(14) njac = iwork(15) nlu = iwork(19) elseif (solver.eq.'VODE') then nsteps = iwork(11)+iwork(21)+iwork(22) naccpt = iwork(11) nfe = iwork(12) njac = iwork(13) nlu = iwork(19) else write(*,*) 'ERROR: report.f: does not support ', + solver, ' (yet)' stop endif ConfigureSolver: not all solvers provide all integration characteristics write(*,*) write(*,'('//LF//'a)') + 'Integration characteristics:' write(*,*) write(*,'('//LF//'3x,a,t35,i8)') + 'number of integration steps', nsteps write(*,'('//LF//'3x,a,t35,i8)') + 'number of accepted steps', naccpt write(*,'('//LF//'3x,a,t35,i8)') + 'number of f evaluations', nfe write(*,'('//LF//'3x,a,t35,i8)') + 'number of Jacobian evaluations', njac if (solver.ne.'DASSL') then write(*,'('//LF//'3x,a,t35,i8)') + 'number of LU decompositions', nlu endif write(*,*) write(*,'('//LF//'a,t38,f8.4,a)') + 'CPU-time used:', cputim, ' sec' if (batch()) then write(*,*) write(*,'(a,'// + '3(a,'//FS//'),'// + '4(e12.3,'//FS//'),'// + 'f8.2,'//FS//','// + '5(i8,'//FS//'),'// + 'f10.4,'//FS//','// + 'f8.2,'//FS//')') + 'TestRun: ', + driver, problm, solver, + rtol, atol, h0, solmax, + scd, + nsteps, naccpt, nfe, njac, nlu, + cputim, + mescd endif return end logical function batch() batch = .false. 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