program integra
c 
c     Author:         Alfio Borzi
c                         Institut für Mathematik
c                         Karl-Franzens-Univrsität Graz
c                         Heinrichstr. 36, A-8010 Graz
c                         Austria
c                         E-mail: alfio.borzi@kfunigraz.ac.at
c
c     Date:            24-01-2000
c     Project:         Multigrid for integral equations 
c     Language:     FORTRAN 77
c     Constraints:   None - Public domain software         
c-----------------------------------------------------------------------------
c     For the presentation of this code see:
c (1)    Alfio Borzi' and Anni Koubek, 
c         Computer Physics Communications 75 (1993) 118-126.
c     For the analysis of the convergence properties of this code see:
c (2)    Alfio Borzi' and Anni Koubek,
c         On a Multi-Grid Algorithm for the TBA Equations
c         In P.W. Hemker and P. Wesseling (Eds.), Multigrid Methods IV,
c             International Series on Numerical Mathematics, Vol. 116,
c             Birkhäuser Verlag, Basel, 1994.
c-----------------------------------------------------------------------------
c     The data structure of this code is as given in:
c     A. Brandt,  Multi-Level Adaptive Solutions to Boundary-Value Problems
c     Mathematics of Computation 31 (1977) 333-390.
c     This code is an implementation of the Nonlinear Multi-Grid Method
c     of the Second Kind (FAS scheme) described in:
c     W. Hackbusch
c     Multi-Grid Methods and Applications, Springer-Verlag, Berlin, 1985. 
c-----------------------------------------------------------------------------
c  solution of a system of NL integral equations with a multi-grid algorithm 
c------------------------multi-grid parameters----------------------------
c  non linear nested iteration is employed         
c  nx0 - number of intervals in the coarsest level k=1
c         if your solution is smooth and the domain smaller than -10 to 10
c         use nx0 = 4, otherwise increase
c  h0  - mesh size of the coarsest level
c         calculated automatically
c  x0  - inf extremum of the interval of integration (x0=-bb)
c         bb is calculated in BOUND
c  u   - initial approximated solution to the equations
c         is composed from known term (TN) and the solution for r -> 0;
c         if this is not adapted for your problem, change it.
c  f_ij - kernel of the integral equation
c                       / bb                  
c   ui(y)=tni(y)+sum_j( |  fij(y,x,u(x))dx) , i=0,..,nsys, j=0,..,nsys
c                       /-bb                   
c     provide the kernel directly in the function F 
c  m   - number of levels     
c  ni0 - number of relaxation in the coarsest level
c  ni1 - number of relaxation before transfer to a coarser grid
c  ni2 - number of relaxation before transfer to a finer grid
c  ng  - if ng=1 v-cycle , if ng=2 w-cycle
c  ngi - number of maximally allowed Multi-Grid cycles 
c  lin - initial level where the approximation is given
c  nrel- if nrel=1 only relaxation done if nrel=0 then nmgm
c  irel- number of maximally allowed relaxation (when nrel=1)
c  nsys- number of equations 
c  the solution, at each level, is stored in q(nsys,qdim)  
c                  m
c  qdim= 2 ( nx0 (2 - 1) + m ) 
c  zero=specifies the value of the residual norm to be reached
c--------------------output--------------------------------------------
c  the output is in OUTPUT.DAT
c  and when iwrite=0 we store
c  the residuals at each computing step in RES.DAT,
c  the CPU time used versus logarithm of the residual norm
c  in TIME.DAT (the subroutine time gives the CPU time used, be aware 
c  to put in the right intrinsic function of your machine) and the 
c  solutions in SOL.DAT. 
c----------------------------------------------------------------------
      implicit real *8 (a-h,o-z)
      common/mass/rm(100)
      common/par/nsys
      common/write/iwrite
      common/coef/pi
      common/zero/zero
      common/rr/r
c     
      external u,tn,f
c     
      qdim=100000
c     
      open(unit=22,status='old',file=
     *     'integra.dat')
      open(unit=23,status='old',file=
     *     'coeff.dat')
      open(unit=24,status='unknown',file=
     *     'res.dat')
      open(unit=25,status='unknown',file=
     *     'output.dat')
      open(unit=26,status='unknown',file=
     *     'time.dat')
      open(unit=27,status='unknown',file=
     *     'sol.dat')
c     
      pi=4.0d0*datan(1.0D0)
c     
      read(22,*) i1,i2
      read(22,*) zero,hx
      read(22,*) nrel,iwrite
      read(22,*) max,step,r0
      close(22)

      do 3 i=1,i2
    3 read(23,*) rm(i)
      close(23)
c
      nsys=i2
      nx0=4
      ni0=4
      ni1=1
      ni2=1
      ng=1
      ngi=3
      irel=50
c     
      l=1
      r=r0
      if(iwrite.eq.0) write(24,106) r
      call bound(bb,f)
      if(iwrite.eq.0) write(24,107) bb
      x0=-bb
      h0=2.0d0*bb/dfloat(nx0)
      xx=dlog(h0/hx)/dlog(2.0d0)
      m=1+nint(xx)
      lin=m-1
      if(iwrite.eq.0) then
      write(24,100) nsys,lin,m
      write(24,101) nx0,h0,hx,zero
      write(24,102) ni0,ni1,ni2
      write(24,103) ng,ngi,nrel,irel
      write(24,104) max,step,r0
      endif
c     
      call nlni(qdim,nx0,nrel,irel,lin,h0,x0,m,ni0,ni1,ni2,ng,ngi,
     *  u,tn,f)
c     
      if(iwrite.eq.0) call printsol(m)     
      call time(second)
      write(25,108) second
c     
  100 format(2x,'nsys=',i3,3x,'lin=',i3,5x,'m=',i3)
  101 format(2x,' nx0=',i3,3x,' h0=',e10.4,2x,'hx=',e10.4,
     *       2x,'zero=',e10.4)
  102 format(2x,' ni0=',i3,3x,'ni1=',i3,2x,' ni2=',i3)
  103 format(2x,'  ng=',i3,3x,'ngi=',i3,2x,'nrel=',i3,2x,'irel=',i3)
  104 format(2x,' max=',i3,2x,'step=',e10.4,4x,'r0=',e10.4)
  106 format(2x,'   r=',e10.4)
  107 format(2x,'  Bb=',e10.4)
  108 format(2x,' total cpu time (secs) ',e10.3)
      close(24)
      close(25)
      close(26)
      close(27)
      stop
      end
c----------------------------------------------------------------------

      subroutine bound(bb,f)
c  calculates the boundary (-bb,bb)
      implicit real *8 (a-h,o-z)
      common/mass/rm(100)
      common/par/nsys
      common/rr/r
      common/coef/pi

      zero=5.0d-16
      x=0.0d0
      y=0.0d0

      bb=0.0d0
      do 20 ne=1,nsys
      do 20 nv=1,nsys
      b=0.0d0
   30 b=b+0.50d0
      z=r*rm(nv)*dcosh(b)
      fb=f(ne,nv,x,y,z)
      if(dabs(fb).lt.zero) goto 21
      goto 30
   21 if(b.gt.bb) bb=b
   20 continue
      return
      end
c----------------------------------------------------------------------

      subroutine crsres(k,krhs,f)
c computes the r.h.s. in the coarser level k
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1),irhs(1)
      common/par/nsys
      call key(k,ist,jj,h,x0)
      call key(krhs,irhs,jj,h,x0)
      io=ist(1)
      ir=irhs(1)
      do 2 ns=1,nsys
      do 1 j=1,jj
      call kappau(k,f,ns,j,sum)
    1 q(ns,ir+j)=q(ns,ir+j)+q(ns,io+j)-sum
    2 continue
      return
      end
c----------------------------------------------------------------------

      real *8 function f(ne,nv,y,x,z)
c  kernel of the integral equations
c  phi(x-y)*log(1+e^(-z))/(2 pi)
      implicit real *8 (a-h,o-z)
      common/rr/r
      common/par/nsys
      common/coef/pi
c
c here compute phi_{ab}(x-y)
c we set phi=-1
c
      phi=-1.
      f=-phi*dlog(1.0d0+dexp(-z))/(pi*2.0d0)
      return
      end
c----------------------------------------------------------------------

      real *8 function tn(n,x)
c  known term of the integral equations 
c  tn=r Ma cosh(x)
      implicit real *8 (a-h,o-z)
      common/mass/rm(100)
      common/coef/pi
      common/rr/r
      tn=r*rm(n)*dcosh(x)
      return
      end
c----------------------------------------------------------------------

      real *8 function u(n,x)
c  initial function for the integral equations
c  u=tn
      implicit real *8 (a-h,o-z)
      common/mass/rm(100)
      common/eps0/eps0(10)
      common/coef/pi
      common/rr/r
      u=r*rm(n)*dcosh(x)
      return
      end
c----------------------------------------------------------------------

      subroutine grdfn(n,jmax,hh,x00)
c  define a jmax array vn
      implicit real *8 (a-h,o-z)
      common/grd/nst(200),jmx(200),h(200),x0
      common/iiqq/iq
      nst(n)=iq
      jmx(n)=jmax
      h(n)=hh
      x0=x00
      iq=iq+jmax
      return
      end
c----------------------------------------------------------------------

      subroutine intadd(k1,k2)
c makes the Coarse Grid correction to the fine-grid solution
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1),irhs(1)
      common/par/nsys
      dimension s(25000)
c     s(qdim/4)
      call key(k1,ist,jjf,hf,x0)
      call key(k2,irhs,jjc,hc,x0)
      ir=irhs(1)
      io=ist(1)
      do 1 ns=1,nsys
      do 2 jc=1,jjc
      jfc=2*jc-1
      s(jc)=q(ns,ir+jc)-q(ns,io+jfc)
    2 continue
      q(ns,io+1)=q(ns,io+1)+s(1)
      q(ns,io+2)=q(ns,io+2)+(5.d0*s(1)+15.d0*s(2)
     +           -5.d0*s(3)+s(4))/16.d0
      jjc2=jjc-2
      do 10 jc=2,jjc2
      jf=2*jc-1
      q(ns,io+jf)=q(ns,io+jf)+s(jc)
      q(ns,io+jf+1)=q(ns,io+jf+1)+(-s(jc-1)+9.d0*s(jc)
     +               +9.d0*s(jc+1)-s(jc+2))/16.d0
   10 continue
      q(ns,io+jjf-2)=q(ns,io+jjf-2)+s(jjc-1)
      q(ns,io+jjf-1)=q(ns,io+jjf-1)+(s(jjc-3)-5.*s(jjc-2)
     +                +15.d0*s(jjc-1)+5.d0*s(jjc))/16.d0
      q(ns,io+jjf)=q(ns,io+jjf)+s(jjc)
    1 continue
      return
      end
c----------------------------------------------------------------------

      subroutine intrp3(kc,kf)
c interpolates the initial solution given at level kc to be a solution
c at level kf
      implicit real *8 (a-h,o-z)
      common q(10,100000),iuf(1),iuc(1)
      common/par/nsys
      call key(kc,iuc,jjc,hc,x0)
      call key(kf,iuf,jjf,hf,x0)
      ic0=iuc(1)
      if0=iuf(1)
      do 1 ns=1,nsys
      q(ns,if0+1)=q(ns,ic0+1)
      q(ns,if0+2)=(5.d0*q(ns,ic0+1)+15.d0*q(ns,ic0+2)
     +  -5.d0*q(ns,ic0+3)+q(ns,ic0+4))/16.d0
      jjc2 = jjc-2
      do 10 jc=2,jjc2
      jf = 2*jc-1
      q(ns,if0+jf) = q(ns,ic0+jc)
      q(ns,if0+jf+1) = (-q(ns,ic0+jc-1)+9.d0*q(ns,ic0+jc)
     +               +9.d0*q(ns,ic0+jc+1)-q(ns,ic0+jc+2))/16.d0
   10 continue
      q(ns,if0+jjf-2) = q(ns,ic0+jjc-1)
      q(ns,if0+jjf-1) = (q(ns,ic0+jjc-3)-5.d0*q(ns,ic0+jjc-2)
     +         +15.d0*q(ns,ic0+jjc-1)+5.d0*q(ns,ic0+jjc))/16.d0
      q(ns,if0+jjf) = q(ns,ic0+jjc)
    1 continue
      return
      end
c----------------------------------------------------------------------

      subroutine kappau(k,f,ne,i,sum)
c integration subroutine
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1)
      common/par/nsys
      call key(k,ist,jj,h,x0)
      sum=0.0d0
      y0=x0
      io=ist(1)
      do 2 ns=1,nsys
      y=(i-1)*h+y0
      x=x0
      z=q(ns,io+1)
      sum0=h*0.50d0*f(ne,ns,y,x,z)
      x=(jj-1)*h+x0
      z=q(ns,io+jj)
      sum0=sum0+h*0.50d0*f(ne,ns,y,x,z)
      do 1 l=2,jj-1
      x=(l-1)*h+x0
      z=q(ns,io+l)
      sum0=sum0+h*f(ne,ns,y,x,z)
    1 continue
      sum=sum+sum0
    2 continue
      return
      end
c----------------------------------------------------------------------

      subroutine key(k,ist,jmax,hh,x00)
c  set ist such that vk(j)=q(.,ist(1)+j)
c  and set jmax=jmx(k) and x00=x0
      implicit real *8 (a-h,o-z)
      common/grd/nst(200),jmx(200),h(200),x0
      dimension ist(1)
      jmax=jmx(k)
      is=nst(k)-jmax-1
      is=is+jmax
      ist(1)=is
      x00=x0
      hh=h(k)
      return
      end

c----------------------------------------------------------------------

      subroutine nlni(qdim,nx0,nrel,irel,lin,h0,x0,m,ni0,ni1,ni2,ng,ngi,
     *   u,tn,f)
c  non-linear nested iteration
      implicit real *8 (a-h,o-z)
      dimension err(10)
      external u,tn,f
      common/par/nsys
      common/write/iwrite
      common/zero/zero
      common/finer/l1
      common/iiqq/iq
      iq=1
      do 1 k=1,m
      k2=2**(k-1)
      call grdfn(k,nx0*k2+1,h0/k2,x0)
    1 call grdfn(k+m,nx0*k2+1,h0/k2,x0)
c
      if(iq.gt.qdim) stop'increase qdim '
c          
      if(iwrite.eq.0) write(26,*)' time, Log(residual norm)'
c
      if(nrel.eq.0) goto 10
      l=m
      call putf(l,u)
      call putf(l+m,tn)
      if(iwrite.eq.0) write(24,*)'level=',l
      if(iwrite.eq.0) write(24,100)
      l1=l
      do 3 i=1,irel
      call relax(l,l+m,f,err)
      errmax=err(1)
      do 5 j=1,nsys
    5 if(err(j).gt.errmax) errmax=err(j)
      call time(second)
      residual=dlog10(errmax)
      if(iwrite.eq.0) write(26,*) second, residual
      if(errmax.lt.zero) goto 11 
    3 continue
      goto 11
c     
   10 call putf(lin,u)
      do 2 l=lin,m
      if(iwrite.eq.0) write(24,*)'level=',l
      if(iwrite.eq.0) write(24,100)
      l1=l
      call nmgm(nx0,nrel,h0,l,m,ni0,ni1,ni2,ng,ngi,u,tn,f)
      if(l.eq.m)goto 11
      call intrp3(l,l+1)
    2 continue
  100 format(3x,'e 1',5x,'e 2',5x,'e 3',5x,'e 4',5x,'e 5'
     *,5x,'e 6',5x,'e 7',5x,'e 8',5x,'e 9',5x,'e 10')
   11 return
      end
c----------------------------------------------------------------------

      subroutine nmgm(nx0,nrel,h0,l,m,ni0,ni1,ni2,ng,ngi,
     *   u,tn,f)
c  non-linear multi-grid method
      implicit real *8 (a-h,o-z)
      dimension err(10)
      common/par/nsys
      common/write/iwrite
      common/zero/zero     
      external u,tn,f
      dimension it(20)
c     
      errmax=0.0d0
      k=l
      call putf(k+m,tn)
      it(k)=ngi
   10 if(k.eq.1) goto 30
   20 do 1 i=1,ni1
      call relax(k,k+m,f,err) 
    1 continue
      if(k.eq.l) then
      errmax=err(1)
      do 60 j=1,nsys
   60 if(err(j).gt.errmax) errmax=err(j)
      if(k.eq.m) then
      call time(second)
      residual=dlog10(errmax)
      if(iwrite.eq.0) write(26,*) second, residual
      end if
      end if
      call rescal(k,k-1+m,k+m,f)
      call putu(k,k-1)
      k=k-1
      call crsres(k,k+m,f)
      it(k)=ng
      go to 10
   30 do 2 i=1,ni0
      call relax(1,1+m,f,err)
    2 continue
   40 if(k.eq.l) return
      k=k+1
      call intadd(k,k-1)
      do 3 i=1,ni2
      call relax(k,k+m,f,err)
    3 continue
      it(k)=it(k)-1
      if(k.eq.l) then
      errmax=err(1)
      do 50 j=1,nsys
   50 if(err(j).gt.errmax) errmax=err(j)
      if(k.eq.m) then
      call time(second)
      residual=dlog10(errmax)
      if(iwrite.eq.0) write(26,*) second, residual
      end if
      end if
      if(errmax.lt.zero.and.k.eq.l) it(l)=0
      if(l.lt.m) it(l)=0
      if(it(k).eq.0) goto 40
      go to 20
      end
c----------------------------------------------------------------------

      subroutine printsol(k)
c  prints (for any nstep) the solutions of the integral equations  
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1)
      common/par/nsys
      common/rr/r     
      call key(k,ist,jj,h,x0)
      nstep=2**5
      jstep=1
      if(jj.gt.(nstep+1)) jstep=(jj-1)/nstep
      write(27,100) r
      write(27,101)
      do 1 n=1,nstep+1
      x=x0+(n-1)*jstep*h
      write(27,102) x,(q(ns,ist(1)+(n-1)*jstep+1),ns=1,nsys)
    1 continue
  100 format(3x,'r= ',e10.4)
  101 format(4x,'x',6x,'s 1',5x,'s 2',5x,'s 3',5x,'s 4',5x,'s 5'
     *,5x,'s 6',5x,'s 7',5x,'s 8',5x,'s 9',5x,'s 10')
  102 format(11(1pe8.1))
      return
      end
c----------------------------------------------------------------------

      subroutine putf(k,f)
c  put the function f at level k
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1)
      common/par/nsys
      call key (k,ist,jj,h,x0)
      io=ist(1)
      do 2 ns=1,nsys
      do 1 j=1,jj
      x=(j-1)*h+x0
    1 q(ns,io+j)=f(ns,x)
    2 continue
      return
      end
c----------------------------------------------------------------------

      subroutine putu(kf,kc)
c  fine-to-coarse grid transfer operator
      implicit real *8 (a-h,o-z)
      common q(10,100000),iuf(1),iuc(1)
      common/par/nsys
      call key(kf,iuf,jjf,hf,x0)
      call key(kc,iuc,jjc,hc,x0)
      if=iuf(1)
      ic=iuc(1)
      do 2 ns=1,nsys
      jf=-1
      do 1 jc=1,jjc
      jf=jf+2
    1 q(ns,ic+jc)=q(ns,if+jf)
    2 continue
      return
      end
c----------------------------------------------------------------------

      subroutine relax(k,krhs,f,err)
c  Gauss-Seidel iteration
      implicit real *8 (a-h,o-z)
      dimension err(10)
      common q(10,100000),ist(1),irhs(1)
      common/finer/m
      common/par/nsys
      common/write/iwrite
c
      call key(k,ist,jj,h,x0)
      call key(krhs,irhs,jj,h,x0)
      ir=irhs(1)
      io=ist(1)
      do 3 ns=1,nsys
      err(ns)=0.0d0
      do 1 j=1,jj
      call kappau(k,f,ns,j,sum)
      qsol=sum+q(ns,ir+j)
c     compute the dynamic residuals
      er=q(ns,io+j)-qsol
      err(ns)=err(ns)+er*er
      q(ns,io+j)=qsol
    1 continue
      err(ns)=dsqrt(err(ns))
    3 continue
      if(k.eq.m.and.iwrite.eq.0) write(24,101)(err(j),j=1,nsys)
  101 format(10(1pe8.1))
      return
      end
c----------------------------------------------------------------------

      subroutine rescal(k,krhc,krhf,f)
c  computes the residual at level k and transfers it on k-1
      implicit real *8 (a-h,o-z)
      common q(10,100000),ist(1),irhc(1),irhf(1)
      common/par/nsys
      call key(k,ist,jjf,hf,x0)
      call key(krhc,irhc,jjc,hc,x0)
      call key(krhf,irhf,jjf,hf,x0)
      io=ist(1)
      if=irhf(1)
      ic=irhc(1)
      do 2 ns=1,nsys
      do 1 jc=1,jjc
      jf=2*jc-1
      call kappau(k,f,ns,jf,sum)
    1 q(ns,ic+jc)=q(ns,if+jf)+sum-q(ns,io+jf)
    2 continue
      return
      end
c----------------------------------------------------------------------

      subroutine time(second)
c  gives the CPU-time
      implicit real *8 (a-h,o-z)
      second=999.
c for unix? use mclock()/100.d0
c for windows use      CALL DCLOCK@(T1)
      return 
      end
c----------------------------------------------------------------------