c----------------------------------------------------------------------- c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c auxiliary routines for all drivers c c DISCLAIMER: see c http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=34 c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~bvpsolvers/src/auxil/report.f c 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 c write(*,*) 'WARNING: report.f: ', c + 'currently no timing; activate timing by configuring', c + 'the function gettim in the file report.f' c gettim = 0 c c if the compiler supports the intrinsic cpu_time 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 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 external lnblnk integer 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 getout(printsolout) logical printsolout c c get output for batch problem c logical batch external batch if (batch()) printsolout=.FALSE. return end subroutine getinp(driver,problm,solver,fullnm, + tolvec,rtol,atol,bvptype) character*(*) driver, problm, solver, fullnm double precision rtol(*), atol(*) logical tolvec character bvptype*5 double precision :: EPS common /pars/ EPS 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 external batch,lnblnk character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) error = .false. ConfigureProblem: straightforward if (.not.( + problm(1:lnblnk(problm)).eq.'bvpT1' .or. + problm(1:lnblnk(problm)).eq.'bvpT1_ho' .or. + problm(1:lnblnk(problm)).eq.'bvpT2' .or. + problm(1:lnblnk(problm)).eq.'bvpT3' .or. + problm(1:lnblnk(problm)).eq.'bvpT4' .or. + problm(1:lnblnk(problm)).eq.'bvpT4_ho' .or. + problm(1:lnblnk(problm)).eq.'bvpT5' .or. + problm(1:lnblnk(problm)).eq.'bvpT6' .or. + problm(1:lnblnk(problm)).eq.'bvpT7' .or. + problm(1:lnblnk(problm)).eq.'bvpT8' .or. + problm(1:lnblnk(problm)).eq.'bvpT9' .or. + problm(1:lnblnk(problm)).eq.'bvpT10' .or. + problm(1:lnblnk(problm)).eq.'bvpT11' .or. + problm(1:lnblnk(problm)).eq.'bvpT12' .or. + problm(1:lnblnk(problm)).eq.'bvpT13' .or. + problm(1:lnblnk(problm)).eq.'bvpT14' .or. + problm(1:lnblnk(problm)).eq.'bvpT15' .or. + problm(1:lnblnk(problm)).eq.'bvpT16' .or. + problm(1:lnblnk(problm)).eq.'bvpT17' .or. + problm(1:lnblnk(problm)).eq.'bvpT18' .or. + problm(1:lnblnk(problm)).eq.'bvpT19' .or. + problm(1:lnblnk(problm)).eq.'bvpT20' .or. + problm(1:lnblnk(problm)).eq.'bvpT21' .or. + problm(1:lnblnk(problm)).eq.'bvpT22' .or. + problm(1:lnblnk(problm)).eq.'bvpT23' .or. + problm(1:lnblnk(problm)).eq.'bvpT24' .or. + problm(1:lnblnk(problm)).eq.'bvpT25' .or. + problm(1:lnblnk(problm)).eq.'bvpT26' .or. + problm(1:lnblnk(problm)).eq.'bvpT27' .or. + problm(1:lnblnk(problm)).eq.'bvpT28' .or. + problm(1:lnblnk(problm)).eq.'bvpT29' .or. + problm(1:lnblnk(problm)).eq.'bvpT30' .or. + problm(1:lnblnk(problm)).eq.'bvpT31' .or. + problm(1:lnblnk(problm)).eq.'bvpT32' .or. + problm(1:lnblnk(problm)).eq.'bvpT33' .or. + problm(1:lnblnk(problm)).eq.'Kidney_Model' .or. + problm(1:lnblnk(problm)).eq.'Measles'.or. + problm(1:lnblnk(problm)).eq.'Flat_Moon' .or. + problm(1:lnblnk(problm)).eq.'FlatEarth' .or. + problm(1:lnblnk(problm)).eq.'FlatEarthDrag' .or. + problm(1:lnblnk(problm)).eq.'Measles' .or. + problm(1:lnblnk(problm)).eq.'Kishalay' + )) +then write(*,*) 'ERROR: report.f: unknown testset problem: ', problm error = .true. endif ConfigureDriver: straightforward if (.not.( + driver.eq.'twpbvpcd' .or. + driver.eq.'twpbvplcd' .or. + driver.eq.'acdccd' .or. + driver.eq.'colmodRd' .or. + driver.eq.'colsysRd' .or. + driver.eq.'colnewRd' .or. + driver.eq.'tomd' .or. + driver.eq.'bvp_m2d' + )) + then write(*,*)'WARNING: report.f: unknown testset driver: ', driver write(*,*) endif ConfigureSolver: straightforward if (.not.( + solver.eq.'TWPBVPC' .or. + solver.eq.'TWPBVPLC' .or. + solver.eq.'ACDCC' .or. + solver.eq.'COLMODR' .or. + solver.eq.'COLSYSR' .or. + solver.eq.'COLNEWR' .or. + solver.eq.'BVP_M2' .or. + solver.eq.'TOM' + )) +then write(*,*) 'ERROR: report.f: unknown testset solver: ', solver error = .true. endif if (error) stop tolvec = .false. rtol(1) = 0d0 atol(1) = 0d0 if (.not. batch()) then write(*,'('//LF//',a,//,'//LF//',a,1x,a,1x,a,1x,a,/)') + ' Test Set for BVP Solvers (release 0.5)', + ' Solving', + fullnm(:lnblnk(fullnm)), + 'using', + solver(:lnblnk(solver)) ConfigureSolver: read rtol/atol ConfigureProblem: read special tolerance write(*,'('//LF//',a)') 'User input:' write(*,*) write(*,'('//LF//',a'//PROMPT//')') + 'give relative error tolerance: ' read *, rtol(1) c write(*,'('//LF//',a'//PROMPT//')') c + 'give absolute error tolerance: ' c read *, atol(1) atol(1) = rtol(1) ConfigureSingularPerturbationProblems, read EPS if (bvptype(1:lnblnk(bvptype)) .eq. 'SPBVP') then write(*,'('//LF//',a'//PROMPT//')') + 'give eps, parameter for the singular perturbation problem: ' read *, EPS endif else ConfigureSolver: generally INDAM BARI only (i.e. normally batch() .eq. .false.) ConfigureSolver: read rtol, atol, n0, same as above read *, rtol(1) atol(1) = rtol(1) if (bvptype(1:lnblnk(bvptype)) .eq. 'SPBVP') then read *, EPS 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) c if (batch()) then ConfigureSolver: generally INDAM BARI only (i.e. normally batch() .eq. .false.) ConfigureSolver: print rtol, atol, n0 ; see read above c write(*,'('//LF//',a,e12.3)') c + 'with relative error tolerance: ', rtol(1) c write(*,'('//LF//',a,e12.3)') c + ' absolute error tolerance: ', atol(1) c c if (bvptype(1:lnblnk(bvptype)) .eq. 'SPBVP') then c write(*,'('//LF//',a,e12.3)') c + ' eps: ', EPS c endif c endif ConfigureProblem: special (vector) tolerances c This is now done in the function settolerances in the c definition of each problem return end subroutine getscd(mescd,scd,neqn,x,yref,y,problm,tolvec,atol,rtol, + nmsh,printout,fileout,nindsol,indsol) integer neqn, nmsh, nindsol, indsol(*) double precision mescd, scd, yref(neqn,*), y(neqn,*), x(nmsh) character*(*) problm, fileout logical tolvec, printout double precision rtol(*), atol(*) c c compute the accuracy of the numerical solution c integer i,j double precision aerr, aerrmx, aerrl2, rerr, rerrmx, rerrl2 double precision mixed_rerr, mixed_rerrmx, min_value double precision scda,scdr logical skipi, skipime character*1 skipc, skipcme parameter (np=6) integer numa, numr, inc, indv(np), ij character*3 FF, LF, PROMPT c c the double precision ieee rounding unit c double precision uround parameter (uround=1.D-16) call getfmt(FF,LF,PROMPT) if (printout) then open(UNIT=90,FILE=fileout) endif numa = 0 numr = 0 numme = 0 aerrl2 = 0d0 rerrl2 = 0d0 aerrmx = 0d0 rerrmx = 0d0 mixed_rerrmx = 0d0 min_value = 1.0d-8 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 inc = (nmsh-1)/np indv(1) = 1 do ij=2,np-1 indv(ij)=indv(ij-1)+inc end do indv(np) = nmsh if (printout) then write(*,'('//LF//',a)') 'Numerical solution:' write(*,*) write(*,'('//LF//',a,t41,a)') +' ' , +' scd ', +' solution component ', +' --------------------------- ignore ', +' + ', +' mixed abs rel mix - abs,rel + ', +'---------------------------------------------------------------', +'--------------- ----- ----- ----- ------------- ' endif do 10 i=1,neqn ij = 1 do 20 j=1,nmsh ConfigureProblem: skipping components in scd computation skipi = .false. skipime = .false. if (skipi) then skipc = 'y' else skipc = ' ' endif if (skipime) then skipcme = 'y' else skipcme = ' ' endif aerr = abs(yref(i,j)-y(i,j)) if (.not. skipi) then aerrmx = max(aerrmx,aerr) aerrl2 = aerrl2 + aerr**2 numa = numa + 1 endif if (abs(yref(i,j)).ge. min_value) then rerr = abs((yref(i,j)-y(i,j))/(yref(i,j))) 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,j)-y(i,j)) + /((atol(i)/rtol(i))+abs(yref(i,j)))) if (.not. skipime) then mixed_rerrmx = max(mixed_rerrmx,mixed_rerr) numme = numme+1 end if if (printout .and. j .eq. indv(ij)) then ij=ij+1 if (aerr.eq.0d0) then write(*, + '('//LF// ', ''x('',i6,'')= '','// + 'e8.2e2,2x,'//' + ''y('',i3,'','',i6,'')='',1x,'// + 'e24.16e3, t66,a,t70,a)' + ) j,x(j),i,j, y(i,j), skipcme,skipc elseif (abs(yref(i,j)).le. min_value ) then write(*, + '('//LF// ',''x('',i6,'')= '','// + 'e8.2e2,2x,'//' + ''y('',i3,'','',i6,'')='',1x,'// + 'e24.16e3,t56,f10.2,1x,f10.2, t66,a,t70,a)' + ) j,x(j),i,j, y(i,j),-log10(mixed_rerr),-log10(aerr), + skipcme, + skipc else write(*, + '('//LF// ', ''x('',i6,'')= '','// + 'e8.2e2,2x,'//' + ''y('',i3,'','',i6,'')='',1x,'// + 'e24.16e3,t56,f10.2,1x,f10.2,1x,f10.2,t66,a,t70,a)' + ) j,x(j),i,j,y(i,j),-log10(mixed_rerr), + -log10(aerr),-log10(rerr), + skipcme, skipc endif end if 20 continue 10 continue if (printout) then write(*,*) nindsol,indsol(1:nindsol) do j=1,nmsh write(90,*) x(j), y(indsol(1:nindsol),j) end do end if aerrl2 = sqrt(aerrl2) rerrl2 = sqrt(rerrl2) if (aerrmx .eq.0d0) aerrmx = uround if (rerrmx .eq.0d0) rerrmx = uround if (mixed_rerrmx.eq.0d0) mixed_rerrmx = uround scda = -log10(aerrmx) scdr = -log10(rerrmx) mescd = -log10(mixed_rerrmx) if (printout) then write(*,*) write(*,'('//LF//',a,t36,i10,1x,i10,1x,i10)') + 'used components for scd', numme,numa, numr write(*,'('//LF//',a,t36,f10.2,1x,f10.2,1x,f10.2)') + 'scd of Y (maximum norm)', + mescd, scda, scdr c if (rerrmx .gt. 0d0) then c write(*,'('//LF//',a,t36,f10.2,1x,f10.2,1x,f10.2)') c + 'scd of Y (maximum norm)', c + -log10(mixed_rerrmx), -log10(aerrmx), -log10(rerrmx) c else c write(*,'('//LF//',a,t36,f10.2,1x,f10.2)') c + 'scd of Y (maximum norm)', c + -log10(mixed_rerrmx), -log10(aerrmx) c end if c write(*,'('//LF//',a,t36,f10.2,1x,f10.2)') c + 'scd of Y (Euclidian norm)', c + -log10(aerrl2), -log10(rerrl2) write(*,*) write(*,*) end if ConfigureProblem: use relative (default) or absolute error ConfigureProblem: for the scd computation if (printout) then write(*,'('//LF//',a,t36,f10.2)') + 'using mixed error yields mescd', mescd end if c if (rerrmx .gt. 0d0) then scd = scdr if (printout) then write(*,'('//LF//',a,t58,f10.2)') + 'using relative error yields scd', scd end if c else c scd = 0 c end if close(90) return end subroutine printsol(neqn,nmsh,x,y,problm, + fileout,nindsol,indsol,printsolout) integer neqn, nmsh, nindsol, indsol(nindsol) double precision x(nmsh),y(neqn,nmsh) character*(*) problm, fileout logical printsolout integer np parameter (np=6) integer i, j, inc, indv(np), ij logical batch external batch character*3 FF, LF, PROMPT call getfmt(FF,LF,PROMPT) c c print of the computed solution c useful for computing the reference solution c if (.not. batch()) then inc = (nmsh-1)/np indv(1) = 1 do ij=2,np-1 indv(ij)=indv(ij-1)+inc end do indv(np) = nmsh write(*,'('//LF//',a)') 'Numerical solution:' write(*,*) do i=1,neqn do ij=1,np j = indv(ij) write(*, + '('//LF// ',''x('',i6,'') = '','// + 'e7.2e2,2x,'//' + ''y('',i3,'','',i6,'') = '','// + 'e24.16e3)' + ) j,x(j),i,j, y(i,j) end do end do endif if (printsolout) then open(UNIT=90,FILE=fileout) do j=1,nmsh write(90,*) x(j), y(indsol(1:nindsol),j) end do close(90) end if return end subroutine report( + driver, problm, solver, + rtol, atol, + iset, cputim, scd,mescd,fmesh, + kappa,kappa1,kappa2,gamma1,sigma) character*(*) driver, problm, solver integer iset(*) double precision rtol, atol real cputim double precision scd,mescd double precision :: EPS common /pars/ EPS c c print integration characteristics including CPU time c integer nsteps, naccpt, nfe, njac, nlu, nc,nss,maxmesh integer fmesh logical batch external batch c c print conditioning parameters c double precision kappa,kappa1,kappa2,gamma1,sigma 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.'TWPBVPC') then nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) nstep = iset(5) nreset = iset(6) maxmesh = iset(7) elseif (solver.eq.'TWPBVPLC') then nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) nstep = iset(5) nreset = iset(6) maxmesh = iset(7) elseif (solver.eq.'ACDCC') then nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) nstep = iset(5) C.... Nc Counts The Total Number Of Continuation Steps Taken. nc = iset(6) C.... Nss Counts The Number Of Successful Continuation Steps Taken. nss = iset(7) maxmesh = iset(8) elseif (solver.eq.'TOM') then nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) nstep = iset(5) nreset = iset(6) maxmesh = iset(7) elseif (solver.eq.'COLMODR') then nstep = 0 nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) C.... Nc Counts The Total Number Of Continuation Steps Taken. nc = iset(5) C.... Nss Counts The Number Of Successful Continuation Steps Taken. nss = iset(6) C.... Nss Counts The Number Of Successful Continuation Steps Taken. maxmesh = iset(7)+1 elseif (solver.eq.'COLSYSR' .or. solver.eq.'COLNEWR') then nstep = iset(5) nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) elseif (solver.eq.'BVP_M2' ) then nstep = 0 nfunc = iset(1) njac = iset(2) nbound = iset(3) njacbound = iset(4) maxmesh = iset(6) else write(*,*) 'ERROR: report.f: does not support ', + solver, ' (yet)' stop endif if (.not. batch()) then ConfigureSolver: not all solvers provide all integration characteristics write(*,*) write(*,'('//LF//',a)') + 'Integration characteristics:' write(*,*) IF (SOLVER .NE. 'COLMODR' .AND. SOLVER .NE. 'BVP_M2') THEN write(*,'('//LF//',3x,a,t35,i8)') + 'number of integration steps', nstep END IF c write(*,'('//LF//',3x,a,t35,i8)') c + 'number of accepted steps', naccpt write(*,'('//LF//',3x,a,t35,i8)') + 'number of f evaluations', nfunc write(*,'('//LF//',3x,a,t35,i8)') + 'number of bc evaluations', nbound write(*,'('//LF//',3x,a,t35,i8)') + 'number of f Jac evaluations', njac write(*,'('//LF//',3x,a,t35,i8)') + 'number of bc Jac evaluations', njacbound IF (SOLVER .eq. 'ACDCC' .or. SOLVER .eq. 'COLMODR') THEN write(*,'('//LF//',3x,a,t65,i8)') + 'Total Number Of Continuation Steps Taken', nc write(*,'('//LF//',3x,a,t65,i8)') + 'Total Number Of Successful Continuation Steps Taken', nss write(*,'('//LF//',3x,a,t65,e10.3)') + 'Final value of Continuation Parameter', EPS END IF IF (SOLVER .NE. 'COLSYSR' .AND. SOLVER .NE. 'COLNEWR') THEN write(*,'('//LF//',3x,a,t35,i8)') + 'Max mesh used', maxmesh END IF write(*,'('//LF//',3x,a,t35,i8)') + 'Final mesh used', fmesh IF (SOLVER .eq. 'ACDCC' .or. SOLVER .eq. 'TWPBVPC' + .or. SOLVER .eq. 'TWPBVPLC' .or. SOLVER .eq. 'TOM' ) THEN c print conditioning parameters write(*,*) write(*,'('//LF//',a)') + 'Conditioning parameters:' write(*,*) write(*,'('//LF//',3x,a,t35,e10.4)') + 'kappa', kappa write(*,'('//LF//',3x,a,t35,e10.4)') + 'kappa1', kappa1 write(*,'('//LF//',3x,a,t35,e10.4)') + 'kappa2', kappa2 write(*,'('//LF//',3x,a,t35,e10.4)') + 'gamma1', gamma1 write(*,'('//LF//',3x,a,t35,e10.4)') + 'sigma', sigma ELSE IF (SOLVER .eq. 'BVP_M2') THEN c print conditioning parameter write(*,*) write(*,'('//LF//',a)') + 'Conditioning parameter:' write(*,*) write(*,'('//LF//',3x,a,t35,e10.4)') + 'kappa', kappa END IF write(*,*) write(*,'('//LF//',a,t38,f8.4,a)') + 'CPU-time used:', cputim, ' sec' endif if (batch()) then write(*,*) write(*,'(a,'// + '3(a,'//FS//'),'// + '3(e12.3,'//FS//'),'// + 'e12.2,'//FS//','// + '6(i10,'//FS//'),'// + '5(e12.4,'//FS//'),'// + 'e12.4,'//FS//','// + 'f12.2,'//FS//')') + 'TestRun: ', + driver, problm, solver, + rtol, atol,eps, + scd, + nsteps, nfunc, njac, nbound, + fmesh,maxmesh, + kappa,kappa1,kappa2,gamma1,sigma, + 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 integer function lrnblnk(string) character*(*) string c c return index of last non blank character in string c integer i do 10 i=1,len(string) if (string(i:i) .ne. ' ') then lrnblnk = i return endif 10 continue lrnblnk = 1 return end subroutine mfileout(problm,solver,filepath,nindsol,indsol) c utility that generate a MATLAB, SCILAB, R and Python function to print the plots c of the solution character problm*15, solver*8 character filemout*140,fileindout*140, filepyout*140 character filepath*100, filesciout*140, filerout*140 integer nindsol, indsol(*) external lnblnk integer i write(filemout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', solver(1:lnblnk(solver)),'.m' write(filesciout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', solver(1:lnblnk(solver)),'.sci' write(filerout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', solver(1:lnblnk(solver)),'.r' write(filepyout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', solver(1:lnblnk(solver)),'.py' write(fileindout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', solver(1:lnblnk(solver)),'.ind' c Matlab open(UNIT=95,FILE=filemout) write(95,*)'% This file is automatically generated by ' write(95,*)'% the BVPtestset file report.f ' write(95,*) 'nindsol = ',nindsol, ';' write(95,*) 'indsol = [' do i=1,nindsol write(95,*) indsol(i) end do write(95,*) '];' write(95,*)'cd ',filepath(1:lnblnk(filepath)) write(95,*) 'addpath(''', filepath(1:lnblnk(filepath)),''')' write(95,*) 'fl=BVPfplot(''', problm(1:lnblnk(problm)), + ''',''', solver(1:lnblnk(solver)),''',indsol,1,...' write(95,*) + '''',filepath(1:lnblnk(filepath)),''');' close(95) open(unit=95,FILE=fileindout) write(95,'(i3)') nindsol close(95) c Scilab File open(UNIT=95,FILE=filesciout) write(95,*)'// This file is automatically generated by ' write(95,*)'// the BVPtestset file report.f ' write(95,*) 'nindsol = ',nindsol, ';' write(95,*) 'indsol = [' do i=1,nindsol write(95,*) indsol(i) end do write(95,*) '];' write(95,*) 'getf(strcat([''',filepath(1:lnblnk(filepath)),''', + ''BVPfplot.sci'']))' write(95,*) 'fl=BVPfplot(''', problm(1:lnblnk(problm)), + ''',''', solver(1:lnblnk(solver)),''',indsol,1,...' write(95,*) + '''',filepath(1:lnblnk(filepath)),''');' close(95) open(unit=95,FILE=fileindout) write(95,'(i3)') nindsol close(95) c ---- R file open(UNIT=95,FILE=filerout) write(95,*)'# This file is automatically generated by ' write(95,*)'# the BVPtestset file report.f ' write(95,*) 'nindsol = ',nindsol write(95,*) 'indsol = c(' write(95,*) indsol(1) do i=2,nindsol write(95,*) ',',indsol(i) end do write(95,*) ')' write(95,*) '#change the working directory to the directory' write(95,*) '#of the script "BVPfplot.R"' write(95,*) 'PATH = getwd()' WRITE(95,*) 'setwd(PATH)' WRITE(95,*) 'PATH = paste(PATH,"/",sep="")' write(95,*) 'file = paste(PATH, "BVPfplot.R", sep="")' write(95,*) 'source(file)' write(95,*) 'fl=BVPfplot(''', problm(1:lnblnk(problm)), + ''',''', solver(1:lnblnk(solver)),''',indsol,1, PATH)' close(95) open(unit=95,FILE=fileindout) write(95,'(i3)') nindsol close(95) c python file c -----------------------------------solve the indentation problem---mm open(UNIT=95,FILE=filepyout) write(95,'(a)')'#This file is automatically generated by ' write(95,'(a)')'#the BVPtestset file report.f ' write(95,'(a)')'import numpy as np' write(95,'(a)')'import BVPfplot as plotter' write(95,'(a)')'PATH = ""' write(95,'(a,a,a)')'PROBLEM = "', problm, '"' write(95,'(a,a,a)')'SOLVER = "', solver, '"' write(95,'(a)')'indsol = np.array([' write(95,*)indsol(1) do i=2,nindsol write(95,*)',',indsol(i) end do write(95,'(a)')'])' write(95,'(a,a,a,a,a)')'fl=plotter.BVPfplot(''', + problm(1:lnblnk(problm)), + ''',''', solver(1:lnblnk(solver)),''',indsol,1, PATH)' close(95) open(unit=95,FILE=fileindout) write(95,'(i3)') nindsol close(95) c-----------------------------------------------------------mm return end