      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-----------------------------------------------------------------------
