c----------------------------------------------------------------------- c Sample driver for the code BIM c----------------------------------------------------------------------- program BIMDd c------------------------------------------- implicit none integer MMAX,lwork,liwork,ltime,lstat parameter(MMAX=1000, & liwork = 40 + MMAX, & lwork = 24 + MMAX*(70 + 3*MMAX) ) parameter(ltime=100,lstat=29) double precision y0(MMAX), y(MMAX), dy(MMAX), yref(MMAX), t0, & tf,t(0:ltime),h0,h,rtol(MMAX),atol(MMAX) integer neqn,mljac,mujac,mlmas,mumas,ind(MMAX),ndisc,itol logical numjac,nummas,consis,tolvec external fevalr, jevalr, mevalr, solout character problm*8 ,fullnm*40,type*3 double precision work(lwork),rpar(1), mescd, scd integer iwork(liwork),ierr,idid,ipar(1),ijac,imas, iout, & icount(lstat),NSTEPS,NACCEPT,NFAILNEWT,NFAILERR,NITER, & i,j real time_begin, time_end, tick, gettim character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) c----------------------------------------------------------------------- c get the problem dependent parameters c----------------------------------------------------------------------- call prob(fullnm,problm,type, neqn,ndisc,t, & numjac,mljac,mujac,nummas,mlmas,mumas,ind) c----------------------------------------------------------------------- c read input values c----------------------------------------------------------------------- if (problm.eq.'pump') then write(6,*) 'give (pump) error tolerance' read(5,*) rtol(1) else write(6,*) 'give relative error tolerance' read(5,*) rtol(1) write(6,*) 'give absolute error tolerance' read(5,*) atol(1) end if write(6,*) 'give the initial stepsize' read(5,*) h0 iwork(9) = 0 iwork(10) = 0 iwork(11) = 0 if (type.eq.'ODE') then imas = 0 elseif (type.eq.'DAE') then imas = 1 do 30 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 *, 'BIMD: ERROR: ', & 'BIMD can not solve index ', ind(i), ' problems' stop end if 30 continue do 40 i=2,neqn if (ind(i).lt.ind(i-1)) then print *, 'BIMD: ERROR: ', & 'BIMD requires the index 1,2,3 variables ', & 'to appear in this order' stop end if 40 continue else print *, 'BIMD: ERROR: ', & 'unknown problem type', type stop end if if (numjac) then ijac = 0 else ijac = 1 end if do i=1,8 iwork(i) = 0 end do do i=1,15 work(i) = 0d0 end do do i=1,lstat icount(i)=0 end do iout = 0 idid = 0 if (problm.eq.'pump') then c charges are much smaller: do i=1,5 atol(i) = rtol(1)*1d-6 rtol(i) = rtol(1) end do do i=6,neqn atol(i) = rtol(1) rtol(i) = rtol(1) end do itol = 1 tolvec = .true. elseif(problm.eq.'water') then do i = 1,36 atol(i) = rtol(1) rtol(i) = rtol(1) end do c pressures are much larger do i = 37,neqn atol(i) = rtol(1)*1d6 rtol(i) = rtol(1) end do itol = 1 tolvec = .true. else itol = 0 tolvec = .false. end if c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(neqn,t(0),y0,dy,consis) do i=1,neqn y(i)=y0(i) end do c----------------------------------------------------------------------- c call of the subroutine BIM c----------------------------------------------------------------------- tick = gettim() tick = gettim() - tick time_begin = gettim() do j=0,ndisc h = h0 t0 = t(j) tf = t(j+1) call BIMD(neqn,fevalr,t0,tf,y,h,rtol,atol,itol, & jevalr,ijac,mljac,mujac, & mevalr,imas,mlmas,mumas, & solout,iout, & work,lwork,iwork,liwork, & rpar,ipar,idid) if (idid.ne.0) then write(6,*) 'ERROR: returned idid =', idid stop endif do i=1,lstat icount(i) = icount(i) + iwork(11+i) end do end do time_end = gettim() c----------------------------------------------------------------------- c compute the integration characteristic c----------------------------------------------------------------------- NSTEPS = 0 NACCEPT = 0 NFAILERR = 0 DO I=1,5 NSTEPS = NSTEPS + icount(I+9) NACCEPT = NACCEPT + icount(I+14) NFAILERR = NFAILERR + icount(I+19) END DO NFAILNEWT = NSTEPS - NACCEPT - NFAILERR c----------------------------------------------------------------------- c print final solution and the integration characteristic c----------------------------------------------------------------------- if (problm.ne.'pump') then write(6,50) atol(1),rtol(1),h0 50 format(//,' we solved the problem with',//, + ' absolute tolerance = ',d10.4,',',/, + ' relative tolerance = ',d10.4,',',/, + ' and initial stepsize = ',d10.4,'.'/) else write(6,55) rtol(1),h0 55 format(//,' we solved the problem with',//, + ' (pump) tolerance = ',d10.4,',',/, + ' and initial stepsize = ',d10.4,'.'/) end if call solut(neqn,t(ndisc+1),yref) call getscd(mescd,scd,neqn,yref,y,problm,tolvec,atol,rtol) write(*,*) write(*,'('//LF//',a)') + 'Integration characteristics:' write(*,*) write(*,'('//LF//'3x,a,t52,i8)') + 'number of integration steps', nsteps write(*,'('//LF//'3x,a,t52,i8)') + 'number of accepted steps', naccept write(*,'('//LF//'3x,a,t52,i8)') + 'number of refused steps (error test)', nfailerr write(*,'('//LF//'3x,a,t52,i8)') + 'number of refused steps (Newton''s convergence)', + nfailnewt write(*,'('//LF//'3x,a,t52,i8)') + 'number of function evaluations', icount(1) write(*,'('//LF//'3x,a,t52,i8)') + 'number of Jacobian evaluations', icount(2) write(*,'('//LF//'3x,a,t52,i8)') + 'number of LU factorization', icount(3) write(*,*) write(*,'('//LF//'a,t38,f8.4,a)') + 'CPU-time used:', time_end-time_begin -tick , ' sec' stop 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) ierr = 0 call feval(neqn,t,y,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) ierr = 0 call jeval(ldjac,neqn,t,y,y,jac,ierr,rpar,ipar) return end subroutine mevalr(neqn,mas,ldmas,ierr,rpar,ipar) double precision mas,rpar(*) integer neqn,ldmas,ierr,ipar(*) dimension y(neqn),mas(ldmas,neqn) double precision t,y ierr = 0 call meval(ldmas,neqn,t,y,y,mas,ierr,rpar,ipar) return end 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(*) integer i,j 20 format (e23.15) write(30,20) t(k) do i=1,m write(30,20) y(i,k) end do write(30,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----------------------------------------------------------------------- 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 mixed_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 c + .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))+abs(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 c --------------------------------------------------------------- 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 c --------------------------------------------------------------------- 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: currently no timing;s activated' gettim = 0 c c if the compiler supports the intrinsic cpu_time c call cpu_time(gettim) c 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