C----------------------------------------------------------------------- c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolve/ c c generic COLMOD driver 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/bvpTestSet/fortransrc/drivers/colmodRd.f c c c----------------------------------------------------------------------- program colmodRd IMPLICIT NONE INTEGER LWRKFL, LWRKIN, NUDIM, NXXDIM, LISERIES, IPRINT, IND, + NTOL, NLBC, NMSH, NFXPNT, NMAX, + IFLBVP, INDNMS, K, NMINIT, IDUM, + NUGDIM DOUBLE PRECISION ALEFT, ARIGHT, CKAPPA1, GAMMA1, + CKAPPA,CKAPPA2,SIGMA, EPS, EPSIN, EPSMIN PARAMETER(NUDIM=20,NXXDIM=10000) Integer MstarI, Kd, Kdm, Nrec Integer NsizeI, NsizeF,NcolPoint parameter (NcolPoint =4) parameter(MstarI = NUDIM, Kd = NcolPoint*nudim, Kdm =Kd+MstarI) parameter(Nrec=10) PARAMETER(Nsizef = NXXDIM*(8 + 4 * MstarI + (5+Kd) * Kdm + + (2*MstarI-Nrec) * 2*MstarI + Kd) ) PARAMETER(NsizeI= (3 + Kdm)*NXXDIM) PARAMETER(NUGDIM=NUDIM) DOUBLE PRECISION FSPACE(NSIZEF) INTEGER ISPACE(NSIZEI) DOUBLE PRECISION U(NUDIM,NXXDIM), XX(NXXDIM) DOUBLE PRECISION YREF(NUDIM,NXXDIM) DOUBLE PRECISION XGUESS(NXXDIM),YGUESS(NUDIM,NXXDIM) DOUBLE PRECISION DMGUESS(NUDIM,NXXDIM) DOUBLE PRECISION UACC(NUDIM,NXXDIM) DOUBLE PRECISION TOL(NUDIM), FIXPNT(1) DOUBLE PRECISION RPAR(10), h DOUBLE PRECISION ATOL(NUDIM), ZETA(NUDIM) double precision yval0(1:nudim) INTEGER LTOL(NUDIM),IPAR(10) INTEGER MS(NUDIM) INTEGER nmguess INTEGER iset(13), icount(8) DOUBLE PRECISION precis(3) integer n0, mstar LOGICAL LINEAR, GIVEU, GIVMSH,giv_u logical full, ACCMESHFAIL LOGICAL use_c, comp_c, tolvec EXTERNAL colmodR EXTERNAL odef EXTERNAL odejac EXTERNAL gsub EXTERNAL dgsub EXTERNAL fnumjac EXTERNAL bcnumjac EXTERNAL MGUESS c this common need to be defined in order to run COLMOD successfully c give information about some common parameters integer nudimval, i common /nudimcom/ nudimval,nmguess common /inguess/ xguess,yguess,dmguess common /neqnval/ neqn common /mstarval/ mstar character fullnm*25, problm*15, bvptype*5 character driver*9, solver*8 parameter (driver = 'colmodRd', solver='COLMODR') real gettim, timtim, cputim external gettim double precision scd, mescd character fileout*140, filepath*100 character namefile*100 logical printsolout, solref, numjac, numbcjac integer nindsol, indsol(nudim),neqn logical stab_kappa, GIVEPS common /pars/ EPS c----------------------------------------------------------------------- c check the problem definition interface date c----------------------------------------------------------------------- call chkdat(driver,solver,20130531) c----------------------------------------------------------------------- c get the problem dependent parameters c----------------------------------------------------------------------- nudimval=nudim fullnm(1:20)=' ' fileout(1:140)=' ' filepath(1:100)=' ' namefile(1:100)=' ' problm(1:12)=' ' bvptype(1:5) = ' ' rpar(1:10) = 0.0d0 ipar(1:10) = 0 call prob(fullnm,problm,bvptype, + neqn,aleft,aright,nlbc,ms, + numjac,numbcjac,linear,rpar,ipar) mstar=0 do i=1,neqn mstar=mstar+ms(i) end do if ( bvptype(1:lnblnk(bvptype)) .ne. 'BVP' .and. + bvptype(1:lnblnk(bvptype)) .ne. 'SPBVP' )then print *, 'COLMODR ERROR: ', + 'unknown Test Set problem type ', bvptype stop endif if ( bvptype(1:lnblnk(bvptype)) .ne. 'SPBVP' ) then print *, 'COLMODRD: ERROR: ', + 'COLMODR cannot solve ', bvptype, ' problems ' stop endif c----------------------------------------------------------------------- c read the tolerances and initial stepsize c----------------------------------------------------------------------- call getinp(driver,problm,solver,fullnm, + tolvec,tol,atol,bvptype) c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(neqn,aleft,aright,ms,yval0,nmsh,givmsh,giveu,xguess, + yguess,dmguess,nudim) n0 = nmsh giv_u=giveu nmguess = nmsh XX(1:NMSH)=XGUESS(1:NMSH) XX(1) = aleft XX(NMSH) = aright IF ( .not. givmsh) then h = (aright-aleft)/(nmsh-1) xx(1) = aleft do i=2,nmsh xx(i)=xx(i-1)+h end do xx(nmsh)= aright END IF call settolerances(neqn,tol,atol,tolvec,ntol,ltol) call setoutput(neqn,solref,printsolout, + nindsol,indsol) call getout(printsolout) EPSMIN = EPS EPS = RPAR(1) DO ind=1,NsizeF Fspace(ind) = 0d0 ENDDO DO ind=1,NsizeI Ispace(ind) = 0 ENDDO c linear not linear IF (LINEAR) THEN iset(1)=0 else iset(1) = 1 end if cNo. Of Collocation Points Per Subinterval ( = K ) iset(2)= 0 Iset(3) = 5 Iset(4) = Ntol Iset(5) = NsizeF Iset(6) = NsizeI c iprint Iset(7) = 1 IF (giveu .and. givmsh) then Iset(8) = 1 Fspace(1:nmguess)=xguess(1:nmguess) Iset(3) = nmguess-1 Iset(9) = 1 else if (givmsh) then Iset(8) = 1 Fspace(1:nmguess)=xguess(1:nmguess) Iset(3) = nmguess-1 Iset(9) = 0 else Iset(8) = 0 Iset(9) = 0 endif Iset(10) = 0 Iset(11) = 0 Iset(12) = 0 IF (RPAR(1) .gt. 0 .and. RPAR(1) > EPSMIN) then Iset(13) = 1 ELSE ISET(13)=0 endif IF (EPS .EQ. 0.0d0) THEN GIVEPS = .FALSE. ELSE GIVEPS = .TRUE. END IF USE_C = .true. COMP_C = .true. precis(3) = epsilon(1.0d0) precis(2) = huge(1.0d0) precis(1) = tiny(1.0d0) IPRINT = -1 DO i=1,nlbc ZETA(i)=Aleft ENDDO DO i=nlbc+1,mstar ZETA(i)=Aright ENDDO NFXPNT= 0 indnms=0 full=.false. if (printsolout) then call getarg(1,filepath) call getarg(2,namefile) if (lnblnk(namefile) .gt. 0) then write(fileout,'(a,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,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 COLMOD c----------------------------------------------------------------------- timtim = gettim() timtim = gettim() - timtim cputim = gettim() if ( bvptype .eq. 'SPBVP') then if (numjac .and. .not. numbcjac) then call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpnt, Ispace, Fspace, Iflbvp, Eps, Epsmin, + odef, fnumjac, gsub, dgsub, mguess, Rpar, Ipar, + Icount) elseif (numjac .and. numbcjac) then call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpnt, Ispace, Fspace, Iflbvp, Eps, Epsmin, + odef, fnumjac, gsub, bcnumjac, mguess, Rpar, Ipar, + Icount) else call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpnt, Ispace, Fspace, Iflbvp, Eps, Epsmin, + odef, odejac, gsub, dgsub, mguess, Rpar, Ipar, + Icount) end if elseif (bvptype .eq. 'BVP') then print *, 'COLMOD: ERROR: ', + 'COLMOD cannot solve ', bvptype, ' problems ' stop end if if (iflbvp.lt.0) then print *, 'COLMOD: ERROR: ', + 'COLMOD returned IFLBVP = ', IFLBVP stop endif if (iflbvp.ne.1) then print *, 'COLMOD: WARNING: ', + 'COLMOD RETURNED A SOLUTION FOR EPSMIN = ', EPS, + 'INSTEAD OF', EPSMIN stop endif nmsh= Ispace(1)+1 DO I=1,NMSH XX(i) = Fspace(I) ENDDO c +Fspace(Is5), Fspace(Is4), DO I=1,NMSH Call MAppsln (XX(i), U(1:neqn,I), Fspace, Ispace) END DO cputim = gettim() - cputim - timtim c----------------------------------------------------------------------- c print numerical solution in endpoint and integration statistics c----------------------------------------------------------------------- c c printsolout = .true. ckappa=0 ckappa1=0 ckappa2=0 gamma1=0 sigma=0 if (solref) then call solut(neqn,xx(1:nmsh),yref(1:neqn,1:nmsh),nmsh,rpar,ipar) call getscd(mescd,scd,neqn,xx(1:nmsh),yref(1:neqn,1:nmsh), + u(1:neqn,1:nmsh), + problm, tolvec,tol,atol, nmsh,printsolout, + fileout,nindsol,indsol) call report( + driver,problm,solver, + tol,tol, + icount,cputim,scd,mescd,nmsh, + ckappa,ckappa1,ckappa2,gamma1,sigma) else c add the computation of a reference solution in the mesh point CALL APPREFSOL(neqn,nmsh,aleft,aright, + xx,u,uacc,tol,ntol,ltol,zeta, + nlbc,MS,iset,linear, fspace, nsizef,ispace,nsizei, + accmeshfail, + nudim,nxxdim,eps,numjac,numbcjac, rpar, ipar) if (.not. accmeshfail) THEN call getscd(mescd,scd,neqn,xx(1:nmsh),uacc(1:neqn,1:nmsh), + u(1:neqn,1:nmsh), + problm, tolvec,tol,atol, nmsh,printsolout, + fileout,nindsol,indsol) call report( + driver,problm,solver, + tol,tol, + icount,cputim,scd,mescd,nmsh, + ckappa,ckappa1,ckappa2,gamma1,sigma) ELSE call printsol(neqn,nmsh,xx(1:nmsh),u(1:neqn,1:nmsh), + problm,fileout,nindsol,indsol,printsolout) end if end if if (printsolout) then close(90) end if end C======================================================================= C `Test Set for BVP Solvers' solout function TWPBVPLC C======================================================================= C======================================================================= C `Test Set for BVP Solvers' ODE wrappers for TWPBVPLC C======================================================================= c c since in TWPBVPLC 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(x,y,f,eps,rpar,ipar) integer ipar(*) double precision x,y(*),f(*),rpar(*) double precision EPS double precision :: EPSVAL integer :: n common /parms/ EPSVAL common /neqnval/ n EPSVAL = EPS call feval(n,x,y(1:n),f(1:n),rpar,ipar) return end subroutine odejac(x,y,dfdy,n,eps,rpar,ipar) implicit none integer n,ipar(*) double precision x,y(*),dfdy(n,*),rpar(*) double precision EPS double precision :: EPSVAL common /parms/ EPSVAL EPSVAL = EPS call jeval(n,x,y,dfdy,rpar,ipar) return end SUBROUTINE GSUB(I,Z,G,EPS,RPAR,IPAR) implicit none integer, INTENT(IN) :: i integer, INTENT(IN) :: ipar(*) double precision, INTENT(IN) :: Z(*) double precision, INTENT(IN) :: rpar(*) double precision, INTENT(OUT) :: G double precision EPS double precision :: EPSVAL common /parms/ EPSVAL integer :: neqn common /neqnval/ neqn EPSVAL = EPS call bceval(I,NEQN,Z,G,rpar,ipar) return end SUBROUTINE DGSUB(I,Z,DG,EPS, RPAR,IPAR) implicit none integer i,ncomp, ipar(*) double precision Z(*), RPAR(*), DG(*) double precision EPS double precision :: EPSVAL common /parms/ EPSVAL common /neqnval/ ncomp EPSVAL = EPS call dbceval(I,NCOMP,Z,DG, RPAR,IPAR) return end c c---------------------------------------------------------------------- c Numerical Jacobian c SUBROUTINE FNUMJAC(X0, Y0, JF0, EPS, RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: IPAR(*) INTEGER :: mstar, R common /mstarval/ mstar DOUBLE PRECISION, INTENT(IN) :: X0, Y0(*), RPAR(*),EPS common /neqnval/ R DOUBLE PRECISION, INTENT(OUT) :: JF0(1:R,*) EXTERNAL odef DOUBLE PRECISION :: DN(MSTAR), Y0N(MSTAR),YSAFE, DELT DOUBLE PRECISION :: X,FP(MSTAR),EPSPREC INTEGER :: I EPSPREC= epsilon(1.0d0) CALL odef(X0,Y0,FP,EPS,RPAR,IPAR) Y0N(1:MSTAR) = Y0(1:MSTAR) DO I=1,MSTAR YSAFE=Y0(I) DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) Y0N(I)=YSAFE+DELT X = X0 CALL odef(X,Y0N,DN,EPS,RPAR,IPAR) JF0(1:R,I)=(DN(1:R)-FP(1:R))/DELT Y0N(I)=YSAFE END DO RETURN END SUBROUTINE FNUMJAC c c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions SUBROUTINE BCNUMJAC(I,y,dbc,EPS,RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: I DOUBLE PRECISION, INTENT(OUT)::dbc(*) DOUBLE PRECISION, INTENT(IN ):: y(*),EPS INTEGER, INTENT(IN):: ipar(*) DOUBLE PRECISION, INTENT(IN) :: rpar(*) INTEGER :: R common /neqnval/ R INTEGER :: mstar common /mstarval/ mstar DOUBLE PRECISION :: DN, Y0N(MSTAR),YSAFE, DELT, BC, EPSPREC INTEGER :: J EXTERNAL gsub EPSPREC= epsilon(1.0d0) Y0N(1:MSTAR) = Y(1:MSTAR) call gsub(i,y,bc,eps,rpar,ipar) DO J=1,MSTAR YSAFE=Y(J) DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) Y0N(J)=YSAFE+DELT call gsub(i,y0N,DN,eps,rpar,ipar) DBC(J)=(DN-BC)/DELT Y0N(J)=YSAFE END DO RETURN END SUBROUTINE BCNUMJAC c c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions SUBROUTINE APPREFSOL(neqn,nmsh,aleft,aright, + xx,u,uacc,tol,ntol,ltol,zeta, + nlbc,MS,iset,linear, fspace, nsizef,ispace,nsizei, + accmeshfail, + nudim,nxxdim,eps,numjac,numbcjac, rpar, ipar) IMPLICIT NONE INTEGER nsizef, nsizei, NUDIM, NXXDIM, LISERIES, IPRINT, IND, + ispace, NTOL, IPAR, LTOL, NLBC, NMSH, NFXPNT, NMAX, + I, NMSHACC, K, NMINIT, IDUM, IFLBVP_A,iset DOUBLE PRECISION TOL, fspace, RPAR, ALEFT, ARIGHT, XX, U, + ZETA, FIXPNTACC, XXACC, UACC,YACC LOGICAL numjac,numbcjac INTEGER MS(NUDIM) DIMENSION FSPACE(*), ISPACE(*), rpar(*),ipar(*),iset(13) DIMENSION U(NUDIM,NXXDIM), XX(NXXDIM),zeta(*) DIMENSION UACC(NUDIM,NXXDIM), XXACC(NXXDIM) DIMENSION YACC(NUDIM,NXXDIM) DOUBLE PRECISION XGUESS(NXXDIM),YGUESS(NUDIM,NXXDIM) DIMENSION TOL(NUDIM), LTOL(NUDIM), FIXPNTACC(NXXDIM) INTEGER indnms_a,iseries_a,nmguess, neqn INTEGER icount_a(8) DOUBLE PRECISION precis(3) LOGICAL LINEAR, GIVEU, GIVMSH logical useC EXTERNAL COLMODR EXTERNAL odef EXTERNAL odejac EXTERNAL gsub EXTERNAL dgsub EXTERNAL fnumjac EXTERNAL bcnumjac EXTERNAL MGUESS logical full, ACCMESHFAIL LOGICAL use_c, comp_c, PRINTSOL, GIVEPS, giv_u integer ureset DOUBLE PRECISION EPS, EPSMIN, EPSIN c this common need to be defined in order to run TWPBVPC successfully c give information about some parameters common/algprs/ nminit, iprint, idum, use_c, comp_c common/acgu/ giv_u, ureset precis(3) = epsilon(1.0d0) precis(2) = huge(1.0d0) precis(1) = tiny(1.0d0) PRINTSOL = .FALSE. full=.false. use_c = .true. useC=use_c NMSHACC= 2*NMSH-1 IF (NMSHACC .gt. NXXDIM ) THEN ACCMESHFAIL = .true. RETURN ENDIF DO I=2,NMSH-1 FIXPNTACC(2*(I-1)-1)=(XX(I-1)+XX(I))/2.0D0 FIXPNTACC(2*(I-1))=XX(I) END DO FIXPNTACC(2*(NMSH-2)+1)=(XX(NMSH-1)+XX(NMSH))/2.0 NFXPNT= 2*(NMSH-2)+1 Xguess(1) = aleft DO I=2,2*(NMSH-2)+2 Xguess(I) =FIXPNTACC(I-1) END DO Xguess(NMSHACC) = aright XXACC(1:NMSHACC)=XGUESS(1:NMSHACC) XXACC(1) = aleft XXACC(NMSHACC) = aright nmguess=nmshacc givmsh = .TRUE. giveu = .TRUE. giveps = .true. Iset(9) = 2 iset(3) = ispace(1) giv_u=giveu epsmin= eps DO I=1,NTOL TOL(I) = TOL(I)/1.0D2 END DO accmeshfail = .false. if (numjac .and. .not. numbcjac) then call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpntacc, Ispace, Fspace, Iflbvp_a, Eps, Epsmin, + odef, fnumjac, gsub, dgsub, mguess, Rpar, Ipar, + Icount_a) elseif (numjac .and. numbcjac) then call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpntacc, Ispace, Fspace, Iflbvp_a, Eps, Epsmin, + odef, fnumjac, gsub, bcnumjac, mguess, Rpar, Ipar, + Icount_a) else call Colmod(neqn, MS, Aleft, Aright, Zeta, Iset, Ltol, + Tol, Fixpntacc, Ispace, Fspace, Iflbvp_a, Eps, Epsmin, + odef, odejac, gsub, dgsub, mguess, Rpar, Ipar, + Icount_a) end if IF (IFLBVP_A .lt. 0 ) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .gt. NXXDIM ) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .eq. NMSH) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .eq. 0 .or. IFLBVP_A .lt. 0) THEN ACCMESHFAIL = .true. RETURN ENDIF DO I=1,NMSH Call MAppsln (XX(i), UACC(1:neqn,I), Fspace, Ispace) END DO RETURN END Subroutine MGuess(X,Z,Dmval,Eps,RPAR,IPAR) double precision X, Z(*),Dmval(*), EPS,RPAR(*) double precision IPAR(*) integer nudim, nxxdim PARAMETER(NUDIM=20,NXXDIM=10000) c nudim must be equal ti nudima c nxxdim must be equal to nxxdim defined in the main fail integer nudima,nmguess,neqn,mstar common /nudimcom/ nudima,nmguess common /inguess/ xguess,yguess,dmguess common /neqnval/ neqn common /mstarval/ mstar double precision xguess(nxxdim), yguess(nudim,nxxdim) double precision dmguess(nudim,nxxdim) double precision xx(1) xx(1)=x call interp(mstar, 1, xx(1), 1, Z, nudima, * nmguess, xguess, yguess) call interp(neqn, 1, xx(1), 1, Dmval, nudima, * nmguess, xguess, dmguess) Return end c =================================================================================== subroutine interp(ncomp, nmsh, xx, nudim, u,nuold, * nmold, xxold, uold) implicit double precision (a-h, o-z) dimension xx(*), u(nudim,*), xxold(*), uold(nuold,*) * blas: dcopy parameter (zero = 0.0d+0) * interp performs piecewise linear interpolation of the old * solution uold at the nmold old mesh points xxold onto the nmsh * new mesh points xx, producing an interpolated solution u. * Note that no assumption is made that the new mesh has * more points than the old, nor that the new and old mesh * points are related in a specific way (except that their first * and last points are identical). * By construction, xx(1) = xxold(1). Copy the first ncomp * components of uold into those of u. call dcopy(ncomp, uold(1,1), 1, u(1,1), 1) i = 2 do 100 im = 2, nmsh-1 50 continue if (i .gt. nmold) return * Check whether the im-th point in the new mesh lies strictly * to the right of, or to the left of (or exactly on) the * i-th point in the old mesh. if (xx(im) .gt. xxold(i)) then i = i + 1 go to 50 else xdif = xxold(i) - xx(im) if (xdif .eq. zero) then * xx(im) and xxold(i) are identical. call dcopy(ncomp, uold(1,i), 1, u(1,im), 1) i = i + 1 else xint = xxold(i) - xxold(i-1) xrat = xdif/xint do 70 k = 1, ncomp u(k,im) = uold(k,i) + xrat*(uold(k,i-1)-uold(k,i)) 70 continue endif endif 100 continue call dcopy(ncomp, uold(1,nmold), 1, u(1,nmsh), 1) return end c ===================================================================================