PROGRAM proof implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/integers/szero,sone,stwo,sthree,sfour common/map/slambda,sc2,sdelta common/banach/smu,snu,snud common/partition/sra,srb,sepspr1,sepspr2,sepsps common/equivct/sequiv_const common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB common/matrix/bm(mps,mps) c********************************************************************* c Resolution: nfp: resolution for the computation of the residual c ---------- n: resolution for the 1st part of the tangent map c m: resolution for the 2nd part of the tangent map c********************************************************************* c Partitions: sra: left point of p_r and p_s c ---------- srb: right point of p_r and p_s c nvb: nbs of intervals in the rough uniform partition Pr c nfv: nbs of intervals in Pr that shall be divided c nf: factor of division c npr1: number of intervals in the left uniform part of p_r c npr2: number of intervals in the right uniform part of p_r c npr1+npr2: total number of intervals in the partition p_r c sepspr1: mesh of the left uniform part of p_r c sepspr2: mesh of the right uniform part of p_r c mps: nbs of intervals in the partiton p_s c nr: nbs of intervals of Pr in each interval of p_s c sepsps: mesh of the uniform partition p_s c********************************************************************* dimension fp(0:1,0:nfp2),fpi(0:1,0:2*nfp2) dimension fp2(0:1,0:nfp2) c****************** c INITIALIZATION * c****************** c Initializes integers szero=siconst(0) sone=siconst(1) stwo=siconst(2) sthree=siconst(3) sfour=siconst(4) c Initializes maps sc2=sdiff(ssqrt(siconst(5)),stwo) slambda_m=srconst(1.7562035d0) slambda_p=srconst(1.7562048d0) sdelta=srconst(rl(squot(slambda_m,slambda_p))) rball=9d-4 c Initializes spaces smu=srconst(5d-1) snu=srconst(9d-1) snud=srconst(ru(squot(snu,sdelta))) c Initializes partitions sra=srconst(6.5d-2) srb=srconst(11.83d0) sepspr2=squot(sdiff(srb,sra),siconst(nvb)) sepspr1=squot(sepspr2,siconst(nf)) sepsps=squot(sdiff(srb,sra),siconst(mps)) print*,' ' c******************************************************************** print*,'COMPUTING THE MATRIX, ITS NORM; SHOWING INVERTIBILITY:' c******************************************************************** call read_fp(fp,0.d0,2) slambda=slambda_p call compute_constant_terms(fp,smu,snu,snud) call compute_matrix call compute_matrix_norm call compute_norm_of_1mSB print*,'Norm of M: ',snorm_of_M print*,'Norm of B: ',snorm_of_B print*,'Norm of (1-S_d)B: ',snorm_of_1mSB call compute_equiv_const print*,'Equivalence constant: ',sequiv_const print*,' ' c************************************************************** print*,'PROVING EXISTENCE FOR lambda in [lambda-,lambda+]:' c************************************************************** c Computes the residual call read_fp(fp,0.d0,2) slambda=slambda_p print*,'lambda+= ',slambda call compute_residual(fp,sepsp,seps,.TRUE.) print*,'Residual for the line: ',seps print*,'Residual for fp+ : ',sepsp c Computes the norm of the tangent map call read_fp(fp,rball,2) call compute_constant_terms(fp,smu,snu,snud) call compute_norm_of_DM(sq) print*,'Norm of DM_(l+,delta): ',ru(sq) c Verifies the existence of the line of fixed point rulhs=ru(seps) rdrhs=rl(sprod(srconst(rball),sdiff(sone,sq))) if(rulhs.lt.rdrhs)then print*,'Existence of the line proved:' print*,rulhs,'<',rdrhs else print*,'Failed to prove existence:' print*,rulhs,'>',rdrhs endif print*,' ' c******************************* print*,'PROVING cl+ > c1:' c******************************* rg=ru(squot(sepsp,sdiff(sone,sq))) call read_fp(fp,rg,2) call fN(fp,fpi,sclp,.TRUE.) sc1=sdiff(sone,sc2) if(rl(sclp).gt.ru(sc1))then print*,'Ok: c1 < cl+' else print*,'Not ok: c1 > cl+' endif print*,'cl+=',sclp print*,'c1 =',sc1 print*,' ' c******************************* print*,'PROVING cl- < c1:' c******************************* c Computes the residual for lambda^- call read_fp(fp,0.d0,1) slambda=slambda_m print*,'lambda-= ',slambda call compute_residual(fp,sepspp,seps,.FALSE.) print*,'Residual for fp- : ',sepspp c Shows that S_delta(fp^-) is in the ball B_rball(fp^+) call fscale_pl(fp,fp2,nfp2,sdelta) call read_fp(fp,0.d0,2) call fdiff(fp2,fp,nfp2,nfp2,fpi) rn=ru(snorm_pl(fpi,2*nfp2,smu,snu)) if(rn.lt.rball)then print*,'S_delta(fp-) in B_r(fp+): ',rn,'<',rball else print*,'S_delta(fp-) not in B_r(fp+): ',rn,'>',rball stop endif c Computes c_lambda- rg=ru(squot(sepspp,sdiff(sone,sq))) fp2(1,nfp2)=sbound(0.d0,rg) slambda=slambda_p call fN(fp2,fpi,sclm,.TRUE.) sclm=sprod(sdelta,sclm) sc1=sdiff(sone,sc2) if(ru(sclm).lt.rl(sc1))then print*,'Ok: cl- < c1' else print*,'Not ok: cl- > c1' endif print*,'cl-=',sclm print*,'c1 =',sc1 print*,' ' end c----------------------------------------------------------------------- c From SECTION 3.1: Interval Analysis c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sbound(rl,ru) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) if(rl.le.ru)then sbound=dcmplx(rl,ru) else print*,'Not an interval: (',rl,',',ru,')' stop endif end c----------------------------------------------------------------------- REAL*8 FUNCTION rl(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) rl=s end c----------------------------------------------------------------------- REAL*8 FUNCTION ru(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) ru=dimag(s) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION siconst(i) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) siconst=sbound(dble(i),dble(i)) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION srconst(rx) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) srconst=sbound(rx,rx) end c----------------------------------------------------------------------- REAL*8 FUNCTION rup(r) implicit real*8(a-h,o-z) logical first_trip save ru,rd,rsru,rsrdn,first_trip data first_trip/.TRUE./ if(first_trip)then r52=1.d0 do i=1,52 r52=r52/2 enddo ru=1.d0+r52 rd=1.d0-r52 print*,'ru=',ru-1.d0 print*,'rd=',1.d0-rd rsru=1.d0 rsrd=1.d0 do i=1,500 rsru=rsru*2 rsrd=rsrd/2 enddo print*,'rsru=',rsru print*,'rsrd=',rsrd rsrdn=-rsrd first_trip=.FALSE. endif if(r.gt.0.d0)then rup=r*ru if(rup.gt.rsru)then print*,'Not a safe bound:',rup stop endif elseif(r.lt.0.d0)then rup=r*rd if(rup.gt.rsrdn)then print*,'Not a safe bound:',rup stop endif else rup=r endif end c----------------------------------------------------------------------- REAL*8 FUNCTION rdown(r) implicit real*8(a-h,o-z) rdown=-rup(-r) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sneg(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) sneg=sbound(-ru(s),-rl(s)) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sabs(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) r=rl(s) sabs=sbound(max(0.d0,r,-ru(s)),max(-r,ru(s))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sinv(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) if(rl(s).gt.0.d0.or.ru(s).lt.0.d0)then sinv=sbound(rdown(1/ru(s)),rup(1/rl(s))) else print*,'Error in sinv, s=',s stop endif end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION spower2(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) sa=sabs(s) rsl=rl(sa) rsu=ru(sa) rpl=rsl*rsl rpu=rsu*rsu resu=rup(max(rpu,rpl)) resl=rdown(min(rpu,rpl)) spower2=sbound(resl,resu) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION ssqrt(s) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) if(rl(s).ge.0.d0)then ssqrt=sbound(rdown(sqrt(rl(s))),rup(sqrt(ru(s)))) else print*,'Error in ssqrt: argument negative ',s stop endif end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION ssum(s1,s2) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) ssum=sbound(rdown(rl(s1)+rl(s2)),rup(ru(s1)+ru(s2))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sdiff(s1,s2) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) sdiff=ssum(s1,sneg(s2)) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sprod(s1,s2) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) r1l=rl(s1) r1u=ru(s1) r2l=rl(s2) r2u=ru(s2) rp1=r1l*r2l rp2=r1l*r2u rp3=r1u*r2l rp4=r1u*r2u resu=rup(max(rp1,rp2,rp3,rp4)) resl=rdown(min(rp1,rp2,rp3,rp4)) sprod=sbound(resl,resu) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION squot(s1,s2) implicit real*8(a-h,o-r,t-z) implicit complex*16(s) squot=sprod(s1,sinv(s2)) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sexp(sx) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour save sxold,sexpold data sxold/(0.d0,0.d0)/ data sexpold/(1.d0,1.d0)/ if(sx.eq.sxold)then sexp=sexpold return endif sxi=sx n=0 2 if(ru(sabs(sxi)).gt.0.03d0)then sxi=squot(sxi,stwo) n=n+1 goto 2 endif sr=sone srx=sxi do i=1,3 sr=ssum(sr,srx) srx=sprod(sxi,squot(srx,siconst(i+1))) enddo re=ru(sprod(sabs(srx),sinv(sdiff(sone,sabs(sxi))))) sr=ssum(sr,sbound(-re,re)) do i=1,n sr=spower2(sr) enddo sxold=sx sexpold=sr sexp=sexpold end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION selogne(seps,n) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour m=4 sr=szero smp=sone smeps=sneg(seps) do k=1,m smp=sprod(smeps,smp) sr=ssum(sr,squot(smp,siconst(n+k))) enddo serr=sabs(squot(sprod(smeps,smp),siconst(n+m+1))) serr=sbound(-ru(serr),ru(serr)) selogne=ssum(sneg(sr),serr) end c----------------------------------------------------------------------- c Discrete Convolution c----------------------------------------------------------------------- SUBROUTINE fft(sdata,isign,n2) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/integers/szero,sone,stwo,sthree,sfour logical first_trip dimension sdata(2*(2*(n2-2))) dimension sang(nfp/2/2+1,2) save sang,first_trip data first_trip/.TRUE./ c Computes angles nx=nfp/2/2+1 if(first_trip)then sco=sinv(ssqrt(stwo)) ssi=sco sang(1,1)=sone sang(1,2)=szero sang(nx,1)=sco sang(nx,2)=ssi nnx=(nx+1)/2 4 sco=ssqrt(squot(ssum(sone,sco),stwo)) ssi=squot(ssi,sprod(stwo,sco)) sang(nnx,1)=sco sang(nnx,2)=ssi do i=2*nnx-1,nx-1,2*nnx-2 sc=sang(i,1) ss=sang(i,2) sang(i+nnx-1,1)=sdiff(sprod(sc,sco),sprod(ss,ssi)) sang(i+nnx-1,2)=ssum(sprod(sc,ssi),sprod(ss,sco)) enddo nnx=(nnx+1)/2 if(nnx.ge.2) goto 4 first_trip=.FALSE. endif c End computing angles nn=2*(n2-2) j=1 do i=1,2*nn,2 if(j.gt.i)then stempr=sdata(j) stempi=sdata(j+1) sdata(j)=sdata(i) sdata(j+1)=sdata(i+1) sdata(i)=stempr sdata(i+1)=stempi endif m=nn 1 if((m.ge.2).and.(j.gt.m))then j=j-m m=m/2 goto 1 endif j=j+m enddo mmax=2 2 if(2*nn.gt.mmax)then istep=2*mmax ist=8*(nx-1)/mmax do m=1,mmax,2 mm=m if(mm.gt.mmax/2) mm=mmax-m+2 k=1+(mm-1)*ist/2 i=1 if(k.gt.nx) then i=2 k=2*nx-k endif swr=sang(k,i) swi=sang(k,3-i) if(isign.eq.-1) swi=sneg(swi) if(m.gt.mmax/2) swr=sneg(swr) do i=m,2*nn,istep j=i+mmax stempr=sdiff(sprod(swr,sdata(j)),sprod(swi,sdata(j+1))) stempi=ssum(sprod(swr,sdata(j+1)),sprod(swi,sdata(j))) sdata(j)=sdiff(sdata(i),stempr) sdata(j+1)=sdiff(sdata(i+1),stempi) sdata(i)=ssum(sdata(i),stempr) sdata(i+1)=ssum(sdata(i+1),stempi) enddo enddo mmax=istep goto 2 endif end c----------------------------------------------------------------------- SUBROUTINE fastconvolution1(f1,sres) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:nfp2) dimension sdata1(2*(2*nfp)) dimension sres(-1:2*nfp2-1) nn=2*nfp do i=1,nfp sdata1(2*i-1)=f1(1,i) sdata1(2*i)=szero enddo do i=2*nfp+1,2*nn sdata1(i)=szero enddo call fft(sdata1,1,nfp2) do i=2,nfp sdir=sdiff(spower2(sdata1(2*i-1)),spower2(sdata1(2*i))) sdim=sprod(stwo,sprod(sdata1(2*i-1),sdata1(2*i))) sdata1(2*i-1)=sdir sdata1(2*i)=sdim sdata1(2*(nn-i+2)-1)=sdata1(2*i-1) sdata1(2*(nn-i+2))=sneg(sdata1(2*i)) enddo do i=1,nfp+1,nfp sdir=sdiff(spower2(sdata1(2*i-1)),spower2(sdata1(2*i))) sdim=sprod(stwo,sprod(sdata1(2*i-1),sdata1(2*i))) sdata1(2*i-1)=sdir sdata1(2*i)=sdim enddo call fft(sdata1,-1,nfp2) sres(-1)=szero sres(0)=szero sres(1)=szero do i=2,nn sres(i)=squot(sdata1(2*(i-1)-1),siconst(nn)) enddo sres(nn+1)=szero sres(nn+2)=szero sres(nn+3)=szero end c----------------------------------------------------------------------- SUBROUTINE fastconvolution2(f1,f2,sres) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(m=1024,m1=m+1,m2=m+2) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:m2),f2(0:1,0:m2) dimension sdata1(2*(2*m)),sdata2(2*(2*m)) dimension sres(-1:2*m2-1) n=m nn=2*m do i=1,m sdata1(2*i-1)=f1(1,i) sdata1(2*i)=szero sdata2(2*i-1)=f2(1,i) sdata2(2*i)=szero enddo do i=2*m+1,2*nn sdata1(i)=szero sdata2(i)=szero enddo call fft(sdata1,1,m2) call fft(sdata2,1,m2) do i=2,m sdir=sdiff(sprod(sdata1(2*i-1),sdata2(2*i-1)), & sprod(sdata1(2*i),sdata2(2*i))) sdim=ssum(sprod(sdata1(2*i-1),sdata2(2*i)), & sprod(sdata1(2*i),sdata2(2*i-1))) sdata1(2*i-1)=sdir sdata1(2*i)=sdim sdata1(2*(nn-i+2)-1)=sdata1(2*i-1) sdata1(2*(nn-i+2))=sneg(sdata1(2*i)) enddo do i=1,m+1,m sdir=sdiff(sprod(sdata1(2*i-1),sdata2(2*i-1)), & sprod(sdata1(2*i),sdata2(2*i))) sdim=ssum(sprod(sdata1(2*i-1),sdata2(2*i)), & sprod(sdata1(2*i),sdata2(2*i-1))) sdata1(2*i-1)=sdir sdata1(2*i)=sdim enddo call fft(sdata1,-1,m2) sres(-1)=szero sres(0)=szero sres(1)=szero do i=2,nn sres(i)=squot(sdata1(2*(i-1)-1),siconst(nn)) enddo sres(nn+1)=szero sres(nn+2)=szero sres(nn+3)=szero end c----------------------------------------------------------------------- c From SECTION 3.2: Standard Sets of B_(alpha beta) c----------------------------------------------------------------------- SUBROUTINE fzero(f,n,sa,seps) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) f(0,n)=seps do i=0,n-1 f(0,i)=ssum(sa,sprod(siconst(i),seps)) enddo do i=0,n f(1,i)=szero enddo end c----------------------------------------------------------------------- SUBROUTINE get_f_on_i(f,n,i,sil,sir,sfl,sfr) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f(0:1,0:n) if((i.ge.1).and.(i.le.n-1))then sil=f(0,i-1) sir=f(0,i) sfl=f(1,i-1) sfr=f(1,i) else print*,'Error in get_f_on_i:',i,n stop endif end c----------------------------------------------------------------------- SUBROUTINE fmult(f1,f2,n,s) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f1(0:1,0:n),f2(0:1,0:n) do i=0,n-1 f2(0,i)=f1(0,i) f2(1,i)=sprod(s,f1(1,i)) enddo f2(0,n)=f1(0,n) f2(1,n)=sprod(sabs(s),f1(1,n)) end c----------------------------------------------------------------------- c From SECTION 4: Operations Involving Functions c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sw(sx,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) sw=sexp(ssum(squot(sal,sx),sprod(sbe,sx))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION ssup_of_w(sd,su,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) swd=sw(sd,sal,sbe) swu=sw(su,sal,sbe) ssup_of_w=srconst(max(ru(swd),ru(swu))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION ssup_of_x_over_w(sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour sx1=ssqrt(ssum(sone,sprod(sfour,sprod(sal,sbe)))) sx2=ssum(sone,sx1) ssup_of_x_over_w=sprod(squot(sx2,sprod(stwo,sbe)),sexp(sneg(sx1))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION ssup_of_winverse(sd,su,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) sxc=ssqrt(squot(sal,sbe)) if(ru(su).lt.rl(sxc))then ssup_of_winverse=srconst(ru(sinv(sw(su,sal,sbe)))) elseif(rl(sd).gt.ru(sxc))then ssup_of_winverse=srconst(ru(sinv(sw(sd,sal,sbe)))) else ssup_of_winverse=srconst(ru(sinv(sw(sxc,sal,sbe)))) endif end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sint_of_w(sd,su,sal,sbe,m,int) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour logical int srd=szero sru=szero seps=squot(sdiff(su,sd),siconst(m)) do i=1,m s1=ssum(sd,sprod(siconst(i-1),seps)) s2=ssum(sd,sprod(siconst(i),seps)) sru=ssum(sru,squot(ssum(sw(s1,sal,sbe),sw(s2,sal,sbe)),stwo)) srd=ssum(srd,sw(squot(ssum(s1,s2),stwo),sal,sbe)) enddo if(int)then sint_of_w=sprod(seps,sbound(rl(srd),ru(sru))) else sint_of_w=squot(sbound(rl(srd),ru(sru)),siconst(m)) endif end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION snorm_pl(f,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) sn=szero do i=1,n-1 call get_f_on_i(f,n,i,sil,sir,sfl,sfr) sv=sprod(sdiff(sir,sil),squot(ssum(sabs(sfr),sabs(sfl)),stwo)) sn=ssum(sn,sprod(sv,ssup_of_w(sil,sir,sal,sbe))) enddo snorm_pl=sn end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION snorm(f,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f(0:1,0:n) sn=snorm_pl(f,n,sal,sbe) snorm=ssum(sn,f(1,n)) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION smass(f,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) if((f(0,n).eq.szero).or. & ((f(1,0).ne.szero).or.(f(1,n-1).ne.szero)))then print*,'Domain error in smass' stop endif s=szero do i=1,n-2 s=ssum(s,f(1,i)) enddo s=sprod(s,f(0,n)) rc=ru(sexp(sneg(sprod(stwo,ssqrt(sprod(sal,sbe)))))) smass=ssum(s,sprod(sbound(-rc,rc),f(1,n))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sexpectation(f,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) if((f(0,n).eq.szero).or. & ((f(1,0).ne.szero).or.(f(1,n-1).ne.szero)))then print*,'Domain error in sexpectation' stop endif seps=f(0,n) s=szero do i=1,n-2 s=ssum(s,sprod(f(0,i),f(1,i))) enddo s=sprod(s,f(0,n)) rc=ru(ssup_of_x_over_w(sal,sbe)) sexpectation=ssum(s,sprod(sbound(-rc,rc),f(1,n))) end c----------------------------------------------------------------------- SUBROUTINE fadd(f1,f2,n1,n2,f3) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:n1),f2(0:1,0:n2),f3(0:1,0:n1+n2) f3(1,n1+n2)=ssum(f1(1,n1),f2(1,n2)) f3(0,n1+n2)=szero i1=0 i2=0 do i=0,n1+n2-1 if(i1.eq.n1)then f3(0,i)=f2(0,i2) f3(1,i)=f2(1,i2) i2=i2+1 elseif(i2.eq.n2)then f3(0,i)=f1(0,i1) f3(1,i)=f1(1,i1) i1=i1+1 elseif(ru(f1(0,i1)).lt.rl(f2(0,i2)))then f3(0,i)=f1(0,i1) if(i2.eq.0)then sv2=szero else call get_f_on_i(f2,n2,i2,si2l,si2r,sf2l,sf2r) sv2=ssum(sf2l,sprod(sdiff(f1(0,i1),si2l), & squot(sdiff(sf2r,sf2l),sdiff(si2r,si2l)))) endif f3(1,i)=ssum(f1(1,i1),sv2) i1=min(i1+1,n1) elseif(ru(f2(0,i2)).lt.rl(f1(0,i1)))then f3(0,i)=f2(0,i2) if(i1.eq.0)then sv1=szero else call get_f_on_i(f1,n1,i1,si1l,si1r,sf1l,sf1r) sv1=ssum(sf1l,sprod(sdiff(f2(0,i2),si1l), & squot(sdiff(sf1r,sf1l),sdiff(si1r,si1l)))) endif f3(1,i)=ssum(f2(1,i2),sv1) i2=min(i2+1,n2) else print*,'Domain error in fadd: overlapping nodes:',i,i1,i2 print*,f1(0,i1),f2(0,i2) stop endif enddo end c----------------------------------------------------------------------- SUBROUTINE fdiff(f1,f2,n1,n2,f3) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f1(0:1,0:n1),f2(0:1,0:n2),f3(0:1,0:n1+n2) do i=0,n2-1 f2(1,i)=sneg(f2(1,i)) enddo call fadd(f1,f2,n1,n2,f3) do i=0,n2-1 f2(1,i)=sneg(f2(1,i)) enddo end c----------------------------------------------------------------------- SUBROUTINE fscale_pl(f1,f2,n,sl) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:n),f2(0:1,0:n) do i=0,n f2(0,i)=squot(f1(0,i),sl) enddo do i=0,n-1 f2(1,i)=sprod(f1(1,i),sl) enddo f2(1,n)=szero end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION fscale_gen(sg,sl,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour if((ru(sl).gt.4d0).or.(ru(sga).gt.rl(sprod(sbe,sl))))then print*,'Domain error in fscal_gen' stop endif s1=sprod(sal,sdiff(sone,squot(sl,sfour))) s2=sprod(s1,sdiff(sbe,squot(sga,sl))) fscale_gen=sprod(sg,sexp(sneg(sprod(stwo,ssqrt(s2))))) end c----------------------------------------------------------------------- SUBROUTINE fscale(f1,f2,n,sl,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f1(0:1,0:n),f2(0:1,0:n) call fscale_pl(f1,f2,n,sl) f2(1,n)=fscale_gen(f1(1,n),sl,sal,sbe,sga) end c----------------------------------------------------------------------- SUBROUTINE ft(f1,f2,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:n),f2(0:1,0:n) c Computes the piecewise affine term f2(0,n)=szero do i=0,n-1 s=sinv(f1(0,i)) f2(0,n-1-i)=srconst((ru(s)+rl(s))/2) enddo do i=0,n-1 s=sprod(f1(1,i),spower2(f1(0,i))) f2(1,n-1-i)=srconst((ru(s)+rl(s))/2) enddo c Computes the general term sstep=f1(0,n) serr=szero do i=1,n-1 call get_f_on_i(f1,n,i,sil,sir,sfl,sfr) sir2=spower2(sir) sil2=spower2(sil) swsi=ssup_of_w(sil,sir,sal,sbe) st1=sprod(selogne(squot(sstep,sil),0),sabs(sdiff(sfr,sfl))) st2=sprod(sstep,sabs(sdiff(squot(sfr,sil),squot(sfl,sir)))) st3=sabs(sdiff(squot(sfr,sil2),squot(sfl,sir2))) st3=sprod(st3,sprod(sstep,sdiff(sir,squot(sstep,stwo)))) serr=ssum(serr,sprod(swsi,ssum(st1,ssum(st2,st3)))) enddo serr=sprod(serr,squot(sstep,sfour)) f2(1,n)=ssum(f1(1,n),serr) end c----------------------------------------------------------------------- SUBROUTINE cubic_spline_coeff(sres,st0,st1,st2,st3,n) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension st0(0:2*n-1),st1(0:2*n-1),st2(0:2*n-1),st3(0:2*n-1) dimension sres(-1:2*n+1) ssix=siconst(6) do i=0,2*n-1 st0(i)=squot(ssum(sres(i+1),ssum(sprod(sres(i),sfour), & sres(i-1))),ssix) st1(i)=squot(sdiff(sres(i+1),sres(i-1)),stwo) st2(i)=sdiff(squot(ssum(sres(i+1),sres(i-1)),stwo), & sres(i)) st3(i)=ssum(squot(sdiff(sres(i+2),sres(i-1)),ssix), & squot(sdiff(sres(i),sres(i+1)),stwo)) enddo end c----------------------------------------------------------------------- SUBROUTINE fcubic_to_pwlinear(st0,st2,st3,f,n,sga,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension st0(0:2*n-1),st2(0:2*n-1),st3(0:2*n-1) dimension f(0:1,0:n+1) sstep=squot(f(0,n+1),stwo) do i=1,n-1 s=sprod(sstep,st0(2*i)) f(1,i)=srconst((rl(s)+ru(s))/2) enddo sfourthird=squot(sfour,sthree) sthreefourth=squot(sthree,sfour) se=szero do i=0,n-1 call get_f_on_i(f,n+1,i+1,sil,sir,sfl,sfr) sei=sprod(ssup_of_w(sil,sir,sga,sbe), & ssum(sprod(sfourthird,sabs(st2(2*i+1))), & sprod(sthreefourth, & ssum(sabs(st3(2*i)),sabs(st3(2*i+1)))))) se=ssum(se,sei) enddo se=sprod(se,spower2(sstep)) f(1,n+1)=ssum(f(1,n+1),se) end c----------------------------------------------------------------------- SUBROUTINE fconvolute1(f1,f3,sal,sga,sbe,comp_se2,se2) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:nfp2),f3(0:1,0:nfp2) dimension st0(0:2*nfp1-1),st1(0:2*nfp1-1), & st2(0:2*nfp1-1),st3(0:2*nfp1-1) dimension sresint(-1:2*nfp2-1) logical comp_se2 if(f1(0,nfp2).eq.szero)then print*,'Domain error in fconvolute1:',f1(0,nfp2) stop endif c Computes the piecewise affine term call fzero(f3,nfp2,sprod(stwo,f1(0,0)),sprod(stwo,f1(0,nfp2))) call fastconvolution1(f1,sresint) call cubic_spline_coeff(sresint,st0,st1,st2,st3,nfp1) call fcubic_to_pwlinear(st0,st2,st3,f3,nfp1,sga,sbe) c Computes the general term sg1=ssum(sprod(stwo,sprod(f1(1,nfp2),snorm_pl(f1,nfp2,sal,sbe))), & spower2(f1(1,nfp2))) f3(1,nfp2)=ssum(f3(1,nfp2),sg1) c Computes E(N^2(f)) if(comp_se2)then se2=sexp_of_tconv(sg1,f1(0,nfp2),f3(0,0), & st0,st1,st2,st3,nfp1,sal,sbe) endif end c----------------------------------------------------------------------- SUBROUTINE fconvolute2(f1,f2,f3,sal,sbe,se2) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(m=1024,m1=m+1,m2=m+2) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:m2),f2(0:1,0:m2),f3(0:1,0:m2) dimension st0(0:2*m1-1),st1(0:2*m1-1),st2(0:2*m1-1),st3(0:2*m1-1) dimension sresint(-1:2*m2-1) r1d=rl(f1(0,m2)) r1u=ru(f1(0,m2)) if((f1(0,m2).ne.f2(0,m2)).or.((r1d.ne.r1u).or.(r1d.eq.0.d0)))then print*,'Domain error in fconvolute2:',f1(0,m2),f2(0,m2) stop endif c Computes the piecewise affine term call fzero(f3,m2,ssum(f1(0,0),f2(0,0)),sprod(stwo,f1(0,m2))) call fastconvolution2(f1,f2,sresint) call cubic_spline_coeff(sresint,st0,st1,st2,st3,m1) call fcubic_to_pwlinear(st0,st2,st3,f3,m1,sal,sbe) c Computes the general term sg1=ssum(sprod(f1(1,m2),snorm(f2,m2,sal,sbe)), & sprod(f2(1,m2),snorm_pl(f1,m2,sal,sbe))) f3(1,m2)=ssum(f3(1,m2),sg1) c Computes E(N^2(f)) se2=sexp_of_tconv(sg1,f1(0,m2),f3(0,0), & st0,st1,st2,st3,m1,sal,sbe) end c----------------------------------------------------------------------- SUBROUTINE fidentity(f1,f2,n1,n2,ra,rlsupp,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/integers/szero,sone,stwo,sthree,sfour dimension f1(0:1,0:n1),f2(0:1,0:n2) dimension fi(0:1,0:2*nfp2) call fzero(f2,n2,srconst(ra),srconst(rlsupp/(n2-1))) call fadd(f1,f2,n1,n2,fi) i=1 do j=0,n1+n2-1 if((fi(0,j).eq.f2(0,i)).and.(i.lt.n2-1))then s=fi(1,j) f2(1,i)=srconst((ru(s)+rl(s))/2) i=i+1 endif enddo call fdiff(f1,f2,n1,n2,fi) f2(1,n2)=ssum(f1(1,n1),snorm_pl(fi,n1+n2,sal,sbe)) end c----------------------------------------------------------------------- SUBROUTINE rsupport(f,n,ra,rlsupp) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/cutoff/rcutoff common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) sa=srconst(ru(squot(ssum(f(0,0),f(0,1)),stwo))) do i=1,n-2 if(abs(ru(f(1,i))).ge.rcutoff)then ra=rl(squot(ssum(f(0,i-1),f(0,i)),stwo)) goto 1 endif enddo print*,'Should not get here; in rsup1:',i,n-2 stop 1 continue do i=n-2,1,-1 if(abs(ru(f(1,i))).ge.rcutoff)then rb=ru(squot(ssum(f(0,i),f(0,i+1)),stwo)) goto 2 endif enddo print*,'Should not get here; in rsup2:',i stop 2 continue rlsupp=rb-ra end c----------------------------------------------------------------------- c From SECTION 5: The Maps N_lambda c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sexp_of_tconv(sg,sstep,sleft, & st0,st1,st2,st3,n,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension st0(0:2*n-1),st1(0:2*n-1),st2(0:2*n-1),st3(0:2*n-1) c Computes (5.13) sk=szero do i=0,2*n-1 sxm=ssum(sleft,sprod(sstep,siconst(i))) seps=squot(sstep,sxm) sk=ssum(sk,sprod(st0(i),selogne(seps,0))) sk=ssum(sk,sprod(st1(i),selogne(seps,1))) sk=ssum(sk,sprod(st2(i),selogne(seps,2))) sk=ssum(sk,sprod(st3(i),selogne(seps,3))) enddo c Computes (5.11) svx=ssup_of_x_over_w(sbe,sprod(sfour,sal)) serrexp=sprod(sbound(-ru(svx),ru(svx)),sg) sexp_of_tconv=ssum(sprod(sstep,sk),serrexp) end c----------------------------------------------------------------------- SUBROUTINE fN(f,fi,scl,eknown) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/map/slambda,sc2,sdelta common/banach/smu,snu,snud common/cutoff/rcutoff common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:nfp2),fi(0:1,0:2*nfp2) dimension f1(0:1,0:nfp2),f2(0:1,0:nfp2) logical eknown c Computes N^1_lambda(f) and N^2_lambda(f) rcutoff=1.d-9 smu4=squot(smu,sfour) call fscale(f,f1,nfp2,slambda,smu,snu,snud) call fconvolute1(f1,f2,smu4,smu,snud,.FALSE.,se2) call ft(f2,f1,nfp2,smu,snud) call rsupport(f1,nfp2,ra,rlsupp) call fidentity(f1,fi,nfp2,nfp2,ra,rlsupp,snud,smu) call fconvolute1(fi,f1,snud,snud,smu,.TRUE.,se2) call ft(f1,fi,nfp2,snud,smu) c Computes c_lambda(f) if(eknown)then se1=squot(sprod(stwo,smass(f,nfp2,smu,snu)),slambda) else se1=squot(sprod(stwo,sprod(sexpectation(f,nfp2,smu,snu), & smass(f,nfp2,smu,snu))),slambda) endif scl=squot(sdiff(sone,sprod(sc2,se2)),se1) c Computes N_lambda(f) call fmult(f2,f1,nfp2,scl) call fmult(fi,f2,nfp2,sc2) call fadd(f1,f2,nfp2,nfp2,fi) end c----------------------------------------------------------------------- SUBROUTINE compute_residual(f,sr1,sr2,both) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/banach/smu,snu,snud common/map/slambda,sc2,sdelta common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB dimension f(0:1,0:nfp2),fi(0:1,0:2*nfp2),fd(0:1,0:3*nfp2) logical both call fN(f,fi,scl,.FALSE.) call fdiff(f,fi,nfp2,2*nfp2,fd) sr1=sprod(snorm_of_M,snorm(fd,3*nfp2,smu,snud)) if(both)then s1=snorm_of_Skappam1(f,nfp2,smu,snu,sdelta) sr2=ssum(sr1,sprod(snorm_of_M,s1)) endif end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION snorm_of_Skappam1(f,n,sal,sbe,skappa) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n) s1=snorm_pl(f,n,sal,sbe) sbk=srconst(ru(squot(sbe,skappa))) s2=szero do i=1,n-1 ssw=ssup_of_w(f(0,i-1),f(0,i),sal,sbk) si2=sprod(sabs(sdiff(f(1,i),f(1,i-1))),ssum(f(0,i-1),f(0,i))) s2=ssum(s2,sprod(ssw,si2)) enddo s2=squot(s2,stwo) snorm_of_Skappam1=sprod(sdiff(sone,skappa),ssum(s1,s2)) end c----------------------------------------------------------------------- c From SECTION 6: The Tangent Maps DN_lambda c----------------------------------------------------------------------- SUBROUTINE compute_constant_terms(fp,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) common/map/slambda,sc2,sdelta common/cutoff/rcutoff common/integers/szero,sone,stwo,sthree,sfour common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp common/tconstantf/fps(0:1,0:n2),fp1(0:1,0:m2),fp1t(0:1,0:m2) dimension fp(0:1,0:nfp2),f1(0:1,0:nfp2) c Computes M(fp), E(fp), and other usefull quantities. smfp=smass(fp,nfp2,sal,sbe) sefp=sexpectation(fp,nfp2,sal,sbe) snfp=snorm(fp,nfp2,sal,sbe) sgfp=fp(1,nfp2) sndfp=snorm_of_der_pl(fp,nfp2,sal,sbe) c Computes N^1_lambda(fp), and other usefull functions. rcutoff=1.d-9 sal4=squot(sal,sfour) call fscale(fp,f1,nfp2,slambda,sal,sbe,sga) call rsupport(f1,nfp2,ra,rlsupp) call fidentity(f1,fps,nfp2,n2,ra,rlsupp,sal4,sga) call fconvolute1(f1,fp,sal4,sal,sga,.FALSE.,se2) snfp1=snorm(fp,nfp2,sal,sga) call rsupport(fp,nfp2,ra,rlsupp) call fidentity(fp,fp1,nfp2,m2,ra,rlsupp,sal,sga) call ft(fp,f1,nfp2,sal,sga) call rsupport(f1,nfp2,ra,rlsupp) call fidentity(f1,fp1t,nfp2,m2,ra,rlsupp,sga,sal) call rsupport(f1,nfp2,ra,rlsupp) call fidentity(f1,fp,nfp2,nfp2,ra,rlsupp,sga,sal) call fconvolute1(fp,f1,sga,sga,sal,.TRUE.,se2) c Computes c_lambda(fp) se1=squot(sprod(stwo,sprod(sefp,smfp)),slambda) scl=squot(sdiff(sone,sprod(sc2,se2)),se1) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION sdelta1(smr,ser,se2) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour common/map/slambda,sc2,sdelta common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp sf1=ssum(squot(smr,smfp),squot(ser,sefp)) sdelta1=sneg((ssum(sprod(scl,sf1), & squot(sprod(slambda,sprod(sc2,se2)), & sprod(stwo,sprod(smfp,sefp)))))) end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION snorm_of_der_pl(f,n2,sal,sbe) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n2) sn=szero do i=1,n2-1 call get_f_on_i(f,n2,i,sil,sir,sfl,sfr) sv=sabs(sdiff(sfr,sfl)) sn=ssum(sn,sprod(sv,ssup_of_w(sil,sir,sal,sbe))) enddo snorm_of_der_pl=sn end c----------------------------------------------------------------------- COMPLEX*16 FUNCTION swsupint(sal,sbe,int) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/partition/sra,srb,sepspr1,sepspr2,sepsps common/integers/szero,sone,stwo,sthree,sfour logical int rrr=0.d0 do i=1,npr1 sd=ssum(sra,sprod(sepspr1,siconst(i-1))) su=ssum(sra,sprod(sepspr1,siconst(i))) s1=sint_of_w(sd,su,sal,sbe,50,int) s2=ssup_of_winverse(sd,su,sal,sbe) rrr=max(rrr,ru(sprod(s1,s2))) enddo do i=1,npr2 sd=ssum(sra,sprod(sepspr2,siconst(nfv+i-1))) su=ssum(sra,sprod(sepspr2,siconst(nfv+i))) s1=sint_of_w(sd,su,sal,sbe,50,int) s2=ssup_of_winverse(sd,su,sal,sbe) rrr=max(rrr,ru(sprod(s1,s2))) enddo swsupint=sbound(0.d0,rrr) end c----------------------------------------------------------------------- SUBROUTINE fDN_center(sn,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour common/map/slambda,sc2,sdelta common/partition/sra,srb,sepspr1,sepspr2,sepsps common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp c Bound (6.17) sg1=ssum(sprod(squot(swsupint(sal,sbe,.TRUE.),stwo),sndfp),sgfp) scf=sprod(stwo,fscale_gen(sg1,slambda,sprod(sfour,sal),sbe,sga)) c Bound (6.11) sinvw=sexp(sneg(sprod(stwo,ssqrt(sprod(sal,sbe))))) rvc=ru(sprod(squot(sepspr2,stwo),sinvw)) ser=sbound(-rvc,rvc) c Bound (6.12) rvx=ru(ssup_of_x_over_w(sal,sprod(sfour,sga))) se2=sprod(stwo,sprod(sbound(-rvx,rvx),sprod(snfp1,scf))) c Bound (6.10) st1=sprod(ssum(sabs(scl),sprod(stwo,sprod(sc2,snfp1))),scf) st2=sprod(sabs(sdelta1(szero,ser,se2)),snfp1) sn=sbound(0.d0,ru(ssum(st1,st2))) end c----------------------------------------------------------------------- SUBROUTINE fDN_left(sn,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour common/map/slambda,sc2,sdelta common/partition/sra,srb,sepspr1,sepspr2,sepsps common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp if((ru(squot(sga,sbe)).gt.rl(slambda)).or. & (ru(slambda).gt.rl(sfour)))then print*,'Domain error in fDN_left' stop endif sxc1=ssqrt(squot(sal,sbe)) sxc2=squot(ssum(sone,ssqrt(ssum(sone,sprod(sfour, & sprod(sal,sbe))))),sprod(stwo,sbe)) if(ru(sra).gt.min(rl(sxc1),rl(sxc2)))then print*,'Domain error in fDN_left' stop endif c Bound (6.24) sx=sprod(ssqrt(slambda),sdiff(stwo,ssqrt(slambda))) scf=sprod(stwo,sprod(sexp(sneg(sprod(sal,squot(sx,sra)))),snfp)) c Bound (6.19) and (6.20) rvc=ru(sinv(sw(sra,sal,sbe))) smr=sbound(-rvc,rvc) rvc=ru(squot(sra,sw(sra,sal,sbe))) ser=sbound(-rvc,rvc) c Bound (6.12) rvx=ru(ssup_of_x_over_w(sal,sprod(sfour,sga))) se2=sprod(stwo,sprod(sbound(-rvx,rvx),sprod(snfp1,scf))) c Bound (6.10) st1=sprod(ssum(sabs(scl),sprod(stwo,sprod(sc2,snfp1))),scf) st2=sprod(sabs(sdelta1(smr,ser,se2)),snfp1) sn=sbound(0.d0,ru(ssum(st1,st2))) end c----------------------------------------------------------------------- SUBROUTINE fDN_right(sn,sal,sbe,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/integers/szero,sone,stwo,sthree,sfour common/map/slambda,sc2,sdelta common/partition/sra,srb,sepspr1,sepspr2,sepsps common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp sxc1=ssqrt(squot(sal,sdiff(sbe,squot(sga,slambda)))) sxc2=squot(ssum(sone,ssqrt(ssum(sone,sprod(sfour, & sprod(sal,sbe))))),sprod(stwo,sbe)) if(rl(srb).lt.max(ru(sxc1),ru(sxc2)))then print*,'Domain error in fDN_right' stop endif c Bound (6.30) sx=sdiff(sbe,squot(sga,slambda)) sexp1=sneg(sprod(stwo,ssqrt(sprod(sal,sx)))) sexp2=sprod(sal,squot(sdiff(slambda,sone),srb)) sexp3=sneg(sprod(srb,sx)) scf1=sprod(stwo,sprod(sexp(ssum(sexp1,ssum(sexp2,sexp3))),snfp)) c Bound (6.31) sexp1=sneg(squot(sprod(sga,srb),slambda)) scf2=sprod(stwo,sprod(sexp(sexp1),sprod(snfp1,scf1))) c Bound (6.32) rvx=ru(ssup_of_x_over_w(sal,sga)) se2=sprod(sbound(-rvx,rvx),scf2) c Bound (6.33) rvc=ru(sinv(sw(srb,sal,sbe))) smr=sbound(-rvc,rvc) rvc=ru(squot(srb,sw(srb,sal,sbe))) ser=sbound(-rvc,rvc) c Bound (6.25) st1=sprod(sabs(scl),scf1) st2=sprod(sc2,scf2) st3=sprod(sabs(sdelta1(smr,ser,se2)),snfp1) sn=sbound(0.d0,ru(ssum(st1,ssum(st2,st3)))) end c----------------------------------------------------------------------- SUBROUTINE fscale_chi(f1,f2,sl) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) dimension f1(3),f2(3) f2(1)=squot(f1(1),sl) f2(2)=squot(f1(2),sl) f2(3)=sprod(f1(3),sl) end c----------------------------------------------------------------------- SUBROUTINE fconv_chi(fb,fp,fi,szeta,sga,seta) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) common/integers/szero,sone,stwo,sthree,sfour dimension fb(3) dimension fp(0:1,0:n2),fi(0:1,0:n2+1) dimension fval(1:n1) save sprevious_delta,fval data sprevious_delta/(0.d0,0.d0)/ c Computes the piecewise affine term if(sprevious_delta.ne.fb(2))then sprevious_delta=fb(2) if(ru(fb(2)).gt.rl(fp(0,n2)))then print*,'Domain error in fconv_chi: delta>eps' stop endif do i=1,n1 sd=squot(sdiff(fp(1,i),fp(1,i-1)),fp(0,n2)) fval(i)=sprod(fb(2),sdiff(fp(1,i),sprod(fb(2),squot(sd,stwo)))) enddo endif call fzero(fi,n2+1,ssum(fb(1),fp(0,0)),fp(0,n2)) do i=1,n1 fi(1,i)=sprod(fb(3),fval(i)) enddo c Computes the general term se=szero do i=1,n sdd=sabs(ssum(fp(1,i+1),sdiff(fp(1,i-1),sprod(fp(1,i),stwo)))) se=ssum(se,sprod(sdd,ssup_of_w(fi(0,i),fi(0,i+1),sga,seta))) enddo se=ssum(se,sprod(sabs(fp(1,1)), & ssup_of_w(fi(0,0),fi(0,1),sga,seta))) se=ssum(se,sprod(sabs(fp(1,n)), & ssup_of_w(fi(0,n1),fi(0,n2),sga,seta))) se1=sprod(spower2(fb(2)),sprod(sdiff(sinv(sfour),squot(fb(2), & sprod(siconst(6),fp(0,n2)))),se)) snfb=sprod(fb(2),ssup_of_w(fb(1),ssum(fb(1),fb(2)),szeta,seta)) se2=sprod(fp(1,n2),snfb) fi(1,n2+1)=sprod(sabs(fb(3)),ssum(se1,se2)) end c----------------------------------------------------------------------- SUBROUTINE fDN_chi(fb,fbi,sal,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) common/map/slambda,sc2,sdelta common/cutoff/rcutoff common/integers/szero,sone,stwo,sthree,sfour common/tconstants/scl,snfp,sndfp,sgfp,snfp1,smfp,sefp common/tconstantf/fps(0:1,0:n2),fp1(0:1,0:m2),fp1t(0:1,0:m2) dimension fb(3),fbs(3) dimension fbi(0:1,0:n2+2*m2+1) dimension f1(0:1,0:n2+1),f2(0:1,0:n2+1) dimension f3(0:1,0:m2),f4(0:1,0:m2),f5(0:1,0:m2) dimension f8(0:1,0:n2+m2+1) rcutoff=1.d-5 sal4=squot(sal,sfour) call fscale_chi(fb,fbs,slambda) call fconv_chi(fbs,fps,f1,sal4,sal,sga) call ft(f1,f2,n2+1,sal,sga) call rsupport(fp1t,m2,ra1,rlsupp1) call rsupport(f2,n2+1,ra2,rlsupp2) if(rlsupp1.gt.rlsupp2)then rlsupp=rlsupp1 else rlsupp=rlsupp2 endif call fidentity(fp1t,f3,m2,m2,ra1,rlsupp,sga,sal) call fidentity(f2,f4,n2+1,m2,ra2,rlsupp,sga,sal) call fconvolute2(f3,f4,f5,sga,sal,se2) call ft(f5,f3,m2,sga,sal) smb=sprod(fb(2),fb(3)) seb=sprod(fb(3),sprod(fb(2),ssum(fb(1),squot(fb(2),stwo)))) se2=sprod(sfour,se2) call fmult(f1,f2,n2+1,sprod(stwo,scl)) call fmult(f3,f4,m2,sprod(sfour,sc2)) call fadd(f2,f4,n2+1,m2,f8) call fmult(fp1,f3,m2,sdelta1(smb,seb,se2)) call fadd(f8,f3,n2+m2+1,m2,fbi) end c----------------------------------------------------------------------- c Frome SECTION 7.1 and 7.2: The Equivalence Constant and the Matrix M c----------------------------------------------------------------------- SUBROUTINE compute_equiv_const implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) common/banach/smu,snu,snud common/equivct/sequiv_const common/integers/szero,sone,stwo,sthree,sfour sequiv_const=ssum(sbound(0.d0,1.d0), & sprod(stwo,swsupint(smu,snu,.FALSE.))) end c----------------------------------------------------------------------- SUBROUTINE projection(f,sv) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/partition/sra,srb,sepspr1,sepspr2,sepsps common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n2+2*m2+1) dimension sv(mps) dimension fr(0:1,0:mps+1),fi(0:1,0:n2+2*m2+mps+2) nm=n2+2*m2+1 call fzero(fr,mps+1,sra,sepsps) call fadd(f,fr,nm,mps+1,fi) if((fi(0,0).ne.fr(0,0)).or.(fi(0,nm+mps).ne.fr(0,mps)))then print*,'Domain error in projection' stop endif i=1 sva=szero do j=0,nm+mps-1 sva=ssum(sva,squot(sprod(sdiff(fi(0,j+1),fi(0,j)), & ssum(fi(1,j),fi(1,j+1))),stwo)) if(fi(0,j+1).eq.fr(0,i))then sv(i)=squot(sva,sepsps) i=i+1 sva=szero endif enddo end c----------------------------------------------------------------------- SUBROUTINE compute_matrix implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/integers/szero,sone,stwo,sthree,sfour common/banach/smu,snu,snud common/partition/sra,srb,sepspr1,sepspr2,sepsps common/matrix/bm(mps,mps) dimension fbase(3) dimension fbasei(0:1,0:n2+2*m2+1) dimension sv(mps) dimension am(mps,mps) dimension bmi(mps,mps) c Computes A --> am do nb=1,mps fbase(1)=ssum(sra,sprod(sepspr2,siconst(nb*nr-nr/2))) fbase(2)=sepspr2 fbase(3)=siconst(nr) call fDN_chi(fbase,fbasei,smu,snud) call projection(fbasei,sv) do i=1,mps am(i,nb)=(ru(sv(i))+rl(sv(i)))/2 enddo enddo c Computes 1-A, --> bmi do i=1,mps do k=1,mps bmi(k,i)=-am(k,i) enddo enddo do i=1,mps bmi(i,i)=1.d0+bmi(i,i) enddo c Computes the numerical inverse of 1-A, --> bmi call gaussj(bmi,mps) c Computes B, --> bm do i=1,mps do j=1,mps x=0.d0 do k=1,mps x=x+am(i,k)*bmi(k,j) enddo bm(i,j)=x enddo enddo c Shows that 1+B is invertible call show_invertibility(bm,am) end c----------------------------------------------------------------------- SUBROUTINE show_invertibility(bm,am) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/integers/szero,sone,stwo,sthree,sfour dimension bm(mps,mps) dimension am(mps,mps) c Computes 1-A do i=1,mps do k=1,mps am(k,i)=-am(k,i) enddo enddo do i=1,mps am(i,i)=1.d0+am(i,i) enddo c Computes the norm of (1+B)(1-A)-1 rmax=0.d0 do i=1,mps sx=szero do j=1,mps sv=szero do k=1,mps if(i.eq.k)then sv=ssum(sv,sprod(ssum(sone,srconst(bm(i,k))), & srconst(am(k,j)))) else sv=ssum(sv,sprod(srconst(bm(i,k)),srconst(am(k,j)))) endif enddo if(i.eq.j) sv=sdiff(sv,sone) sx=ssum(sx,sabs(sv)) enddo rmax=max(rmax,ru(sx)) enddo if(rmax.lt.1.d0)then print*,'Invertible:',rmax,'<1' else print*,'Not invertible:',rmax stop endif end c----------------------------------------------------------------------- SUBROUTINE gaussj(a,n) implicit integer(n,m) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) integer i,icol,irow,j,k,l,ll,indxc(mps),indxr(mps),ipiv(mps) real*8 a(n,n) real*8 big,dum,pivinv do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0.d0 do 13 j=1,n if(ipiv(j).ne.1)then do 12 k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big)then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then print*,'Singular matrix in gaussj 1' stop endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.d0)then print*,'Singular matrix in gaussj 2' stop endif pivinv=1.d0/a(icol,icol) a(icol,icol)=1.d0 do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue do 21 ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0.d0 do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue endif 21 continue 22 continue do 24 l=n,1,-1 if(indxr(l).ne.indxc(l))then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return END c----------------------------------------------------------------------- SUBROUTINE compute_matrix_norm implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/banach/smu,snu,snud common/partition/sra,srb,sepspr1,sepspr2,sepsps common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB common/integers/szero,sone,stwo,sthree,sfour common/matrix/bm(mps,mps) rmax=0.d0 do j=1,mps sx=szero do i=1,mps sd=ssum(sra,sprod(sepsps,siconst(i-1))) su=ssum(sra,sprod(sepsps,siconst(i))) sx=ssum(sx,sprod(sabs(srconst(bm(i,j))), & sint_of_w(sd,su,smu,snu,50,.FALSE.))) enddo sd=ssum(sra,sprod(sepsps,siconst(j-1))) su=ssum(sra,sprod(sepsps,siconst(j))) rmax=max(rmax,ru(sprod(sx,ssup_of_winverse(sd,su,smu,snu)))) enddo snorm_of_B=sbound(0.d0,rmax) snorm_of_M=ssum(sbound(0.d0,1.d0),sbound(0.d0,rmax)) end c----------------------------------------------------------------------- c From SECTION 7.3: The Tangent Maps DM_(lambda,kappa) c----------------------------------------------------------------------- COMPLEX*16 FUNCTION snorm_add(f,sv,sal,sga) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/partition/sra,srb,sepspr1,sepspr2,sepsps common/integers/szero,sone,stwo,sthree,sfour dimension f(0:1,0:n2+2*m2+1),sv(mps) dimension fi(0:1,0:n2+2*m2+mps+2),fr(0:1,0:mps+1) nm=n2+2*m2+1 call fzero(fr,mps+1,sra,sepsps) call fadd(f,fr,nm,mps+1,fi) if((fi(0,0).ne.fr(0,0)).or.(fi(0,nm+mps).ne.fr(0,mps)))then print*,'Domain error in snorm_add' stop endif sn=szero i=1 do j=1,nm+mps call get_f_on_i(fi,nm+mps+1,j,sil,sir,sfl,sfr) sfl=sabs(ssum(sfl,sv(i))) sfr=sabs(ssum(sfr,sv(i))) sva=sprod(sdiff(sir,sil),ssum(sfr,sfl)) sn=ssum(sn,sprod(sva,ssup_of_w(sil,sir,sal,sga))) if(fi(0,j).eq.fr(0,i))then i=i+1 endif enddo snorm_add=sbound(0.d0,ru(squot(sn,stwo))) end c----------------------------------------------------------------------- SUBROUTINE linear_app(sv1,sv2) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/integers/szero,sone,stwo,sthree,sfour common/matrix/bm(mps,mps) dimension sv1(mps) dimension sv2(mps) do i=1,mps sx=szero do j=1,mps sx=ssum(sx,sprod(srconst(bm(i,j)),sv1(j))) enddo sv2(i)=sx enddo end c----------------------------------------------------------------------- SUBROUTINE compute_norm_of_1mSB implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/map/slambda,sc2,sdelta common/banach/smu,snu,snud common/partition/sra,srb,sepspr1,sepspr2,sepsps common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB common/integers/szero,sone,stwo,sthree,sfour rrr=0.d0 do i=1,mps sd=ssum(sra,sprod(sepsps,siconst(i-1))) su=ssum(sra,sprod(sepsps,siconst(i))) swi=sint_of_w(sd,su,smu,snu,50,.TRUE.) sdk=squot(sd,sdelta) suk=squot(su,sdelta) s1=sint_of_w(sd,sdk,smu,snu,50,.TRUE.) s2=sint_of_w(su,suk,smu,snu,50,.TRUE.) rrr=max(rrr,ru(squot(ssum(s1,s2),swi))) enddo rrr=ru(ssum(sdiff(sone,sdelta),srconst(rrr))) snorm_of_1mSB=sprod(snorm_of_B,srconst(rrr)) end c----------------------------------------------------------------------- SUBROUTINE init_chi(fb,svb,nb) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) common/banach/smu,snu,snud common/integers/szero,sone,stwo,sthree,sfour common/partition/sra,srb,sepspr1,sepspr2,sepsps dimension fb(3) if(nb.le.npr1)then sd=ssum(sra,sprod(sepspr1,siconst(nb-1))) su=ssum(sra,sprod(sepspr1,siconst(nb))) sv=sinv(sint_of_w(sd,su,smu,snu,50,.TRUE.)) fb(2)=sepspr1 else sd=ssum(sra,sprod(sepspr2,siconst(nb-npr1+nfv-1))) su=ssum(sra,sprod(sepspr2,siconst(nb-npr1+nfv))) sv=sinv(sint_of_w(sd,su,smu,snu,50,.TRUE.)) fb(2)=sepspr2 endif fb(1)=sd fb(3)=sone svb=sbound(0.d0,ru(sv)) end c----------------------------------------------------------------------- SUBROUTINE fDM_chi(nbase,snormi) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(n=512,n1=n+1,n2=n+2) parameter(m=1024,m1=m+1,m2=m+2) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) parameter(nr=10,mps=nvb/nr) common/banach/smu,snu,snud common/map/slambda,sc2,sdelta common/partition/sra,srb,sepspr1,sepspr2,sepsps common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB common/integers/szero,sone,stwo,sthree,sfour dimension fbase(3) dimension fbasei(0:1,0:n2+2*m2+1) dimension sv1(mps),sv2(mps) c Computes DN on basis vector nbase call init_chi(fbase,svb,nbase) call fDN_chi(fbase,fbasei,smu,snud) sDNpl=snorm_pl(fbasei,n2+2*m2+1,smu,snu) sDNg=fbasei(1,n2+2*m2+1) sSkm1DNpl=snorm_of_Skappam1(fbasei,n2+2*m2+1,smu,snu,sdelta) c Computes (B DNpl) - (B chi) call projection(fbasei,sv1) if(nbase.le.npr1)then k=1+(nbase-1)/nf/nr else k=1+(nbase-npr1+nfv-1)/nr endif smb=squot(sprod(fbase(2),fbase(3)),sepsps) sv1(k)=sdiff(sv1(k),smb) call linear_app(sv1,sv2) c Computes the norm st1=sprod(snorm_of_M,sDNg) st2=snorm_add(fbasei,sv2,smu,snud) st3=sprod(snorm_of_B,sSkm1DNpl) st4=sprod(snorm_of_1mSB,ssum(sprod(svb,sDNpl),sone)) snormi=ssum(sprod(svb,ssum(st1,ssum(st2,st3))),st4) end c----------------------------------------------------------------------- SUBROUTINE compute_norm_of_DM(sq) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nvb=5000,nfv=50,nf=2,npr1=nfv*nf,npr2=nvb-nfv) common/banach/smu,snu,snud common/equivct/sequiv_const common/matrixnorm/snorm_of_M,snorm_of_B,snorm_of_1mSB rn=0.d0 call fDN_center(src,smu,snu,snud) call fDN_left(srl,smu,snu,snud) call fDN_right(srr,smu,snu,snud) sf=sprod(snorm_of_M,sequiv_const) rnc=ru(sprod(sf,src)) rnl=ru(sprod(sf,srl)) rnr=ru(sprod(sf,srr)) print*,'From fDN_center:',rnc print*,'From fDN_left :',rnl print*,'From fDN_right :',rnr rn1=max(rnc,max(rnl,rnr)) rn2=0.d0 do nbase=1,npr1+npr2,1 call fDM_chi(nbase,snormi) rn2=max(rn2,ru(snormi)) enddo rn2=ru(sprod(sequiv_const,sbound(0.d0,rn2))) print*,'From fDM_chi :',rn2 sq=sbound(0.d0,max(rn1,rn2)) end c----------------------------------------------------------------------- c Read Input Data Files c----------------------------------------------------------------------- SUBROUTINE read_fp(f,rgen,k) implicit real*8(a-e,g-h,o-r,t-z) implicit complex*16(s,f) parameter(nfp=131072,nfp1=nfp+1,nfp2=nfp+2) dimension f(0:1,0:nfp2) logical fhere if(k.eq.1)then inquire(file='fpoint.lm',EXIST=fhere) if(fhere)then open(unit=17,status='OLD',file='fpoint.lm') else print*,'File fpoint.lm missing' stop endif else inquire(file='fpoint.lp',EXIST=fhere) if(fhere)then open(unit=17,status='OLD',file='fpoint.lp') else print*,'File fpoint.lp missing' stop endif endif read(17,'(1P,E16.10)',ERR=2)ra read(17,'(1P,E16.10)',ERR=2)rb sfpa=srconst(ra) seps=srconst((rb-ra)/nfp1) call fzero(f,nfp2,sfpa,seps) f(1,nfp2)=sbound(0.d0,rgen) do i=1,nfp read(17,'(1P,E16.10)',IOSTAT=ios,ERR=2)rv f(1,i)=srconst(rv) enddo read(17,*,IOSTAT=ios1)rv if((ios.ne.0).or.(ios1.ne.-1))goto 2 goto 1 2 continue if(k.eq.1)then print*,'Improper file fpoint.lm',ios,ios1 else print*,'Improper file fpoint.lp',ios,ios1 endif stop 1 continue close(17) end c-----------------------------------------------------------------------