




















































































































 


































































































































































 






























































































































































































 





















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                




































c
      subroutine mgparm(m,ifd59,id5,id9,idi,idg,imx,jmx)
      implicit REAL*8 (a-h,o-z)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
      parameter(n5=1,n9=2,ni=3,jm=4,i9=5,j9=6,ifd=7,jred=8)
      dimension np2(20,8)
      iq5=1
      iq9=1
      iqi=1
      m=1
      np2(m,1)=jmx
      np2(m,2)=3
   10 if(np2(m,1).le.3) go to 20
      m=m+1
      np2(m,1)=np2(m-1,1)/2+1
      if(np2(m-1,2).eq.2.and.mod(np2(m-1,1),2).eq.1)
     + np2(m,1)=np2(m,1)+1
      np2(m,2)=2
      go to 10
   20 do 30 k=1,m
      np2(m-k+1,jm)=np2(k,1)
   30 np2(m-k+1,jred)=np2(k,2)
      do 40 k=m,1,-1
      ktot=imx*np2(k,jm)
      np2(k,n5)=iq5
      iq5=iq5+ktot
      np2(k,n9)=iq9
      if(k.lt.m.or.ifd59.eq.9) iq9=iq9+ktot
      np2(k,ni)=iqi
   40 if(k.lt.m) iqi=iqi+ktot
      do 50 k=1,m
      np2(k,i9)=imx
      np2(k,j9)=np2(k,jm)
   50 np2(k,ifd)=9
         if(ifd59.eq.5) then
      np2(m,i9)=1
      np2(m,j9)=1
      np2(m,ifd)=5
         endif
      id5=iq5-1
      id9=iq9-1
      idi=iqi-1
      idg=imx
      return
      end
c
c
c
c
c
      subroutine mg9 (idim,ilower,iupper,jdim,jlower,jupper,
     &             cc,cn,cs,cw,ce,cnw,cne,csw,cse,u,rhs,
     &             id5,id9,idi,idg,ifmg,eps,rmax,ier)
      implicit REAL*8 (a-h,o-z)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
      integer id5,id9,idi,idg,ier
      REAL*8 cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
     &     ce(idim,jdim),cnw(idim,jdim),cne(idim,jdim),csw(idim,jdim),
     &     cse(idim,jdim),u(idim,jdim),rhs(idim,jdim)
      REAL*8 eps,rmax

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      integer ncyc,ifd59
      parameter (ncyc=40,ifd59=9)

      integer np2(20,8)
      integer iskip,ipc,nman
      parameter (iskip=0,ipc=1,nman=0)
      integer irc,irurd,im,jm
      integer linp,lout
      common /io/ linp,lout

      REAL*8 ac(id5),an(id5),as(id5),aw(id5),ae(id5),
     &     anw(id9),ane(id9),asw(id9),ase(id9),
     &     q(id5),f(id5),gam(idg)

      REAL*8 ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)

      REAL*8 tol

      integer iscale,itry


      ier=0

c
c
c
      lout=6
c
      irc=0
      irurd=0
      im=iupper-ilower+3
      jm=jupper-jlower+3
      tol=eps

c
c
c
      do 110 j=jlower,jupper
      do 100 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      ac(n)=cc(i,j)
      an(n)=cn(i,j)
      as(n)=cs(i,j)
      aw(n)=cw(i,j)
      ae(n)=ce(i,j)
      anw(n)=cnw(i,j)
      ane(n)=cne(i,j)
      asw(n)=csw(i,j)
      ase(n)=cse(i,j)
      q(n)=u(i,j)
      f(n)=rhs(i,j)
  100 continue
  110 continue


c
c
c
c
c
      iscale=1
      do 200 j=jlower,jupper
      do 205 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      if (ac(n) .eq. 0.) then
        iscale=0
      endif
  205 continue
  200 continue

c
c
c
      if (iscale.eq.1) then

      do 210 j=jlower,jupper
      do 215 i=ilower,iupper
        n=(j+1-jlower)*im + i+1-(ilower-1)
        f(n)=f(n)/ac(n)
        ase(n)=ase(n)/ac(n)
        asw(n)=asw(n)/ac(n)
        ane(n)=ane(n)/ac(n)
        anw(n)=anw(n)/ac(n)
        ae(n)=ae(n)/ac(n)
        aw(n)=aw(n)/ac(n)
        as(n)=as(n)/ac(n)
        an(n)=an(n)/ac(n)
        ac(n)=1.
  215 continue
  210 continue

      endif

c
c

      itry=0
      
 1122 call umgs2(
     . ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
     . ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
     . ipc,irc,irurd)

      if ((rmax.gt.tol).and.(itry.le.5)) then
        itry=itry+1
        print*,"Retry #",itry
        goto 1122
      endif

      if (rmax.gt.tol) then
        ier = -1
        print*,"Did not converge"
        print*," maximum residual    = ",rmax
        print*," tolerance           = ",tol
      endif

c
c
c
      do 510 j=jlower,jupper
      do 500 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      u(i,j)=q(n)
  500 continue
  510 continue

      return
      end



c
c
c
c
c
      subroutine mg5 (idim,ilower,iupper,jdim,jlower,jupper,
     &             cc,cn,cs,cw,ce,u,rhs,
     &             id5,id9,idi,idg,ifmg,eps,rmax,ier)
      implicit REAL*8 (a-h,o-z)
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
      integer id5,id9,idi,idg,ier
      REAL*8 cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
     &     ce(idim,jdim),u(idim,jdim),rhs(idim,jdim)
      REAL*8 eps,rmax

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      integer ncyc,ifd59
      parameter (ncyc=40,ifd59=5)

c
c
c
c
c
c
c

      integer np2(20,8)
      integer iskip,ipc,nman
      parameter (iskip=0,ipc=1,nman=0)
      integer irc,irurd,im,jm
      integer linp,lout
      common /io/ linp,lout

      REAL*8 ac(id5),an(id5),as(id5),aw(id5),ae(id5),
     &     anw(id9),ane(id9),asw(id9),ase(id9),
     &     q(id5),f(id5),gam(idg)

      REAL*8 ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)

      REAL*8 tol

      integer iscale,itry

c
      ac = 0.
      an = 0.
      as = 0.
      aw = 0.
      ae = 0.
      q = 0.
      f = 0.
      gam = 0.

      ier=0
c
c
c
      lout=6
c
      irc=0
      irurd=0
      im=iupper-ilower+3
      jm=jupper-jlower+3
      tol=eps


c
c
c
      do 110 j=jlower,jupper
      do 100 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      ac(n)=cc(i,j)
      an(n)=cn(i,j)
      as(n)=cs(i,j)
      aw(n)=cw(i,j)
      ae(n)=ce(i,j)
      anw(n)=0.
      ane(n)=0.
      asw(n)=0.
      ase(n)=0.
      q(n)=u(i,j)
      f(n)=rhs(i,j)
  100 continue
  110 continue


c
c
c
c
c
      iscale=1
      do 200 j=jlower,jupper
      do 205 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      if (ac(n) .eq. 0.) then
        iscale=0
      endif
  205 continue
  200 continue

c
c
c
      if (iscale.eq.1) then

      do 210 j=jlower,jupper
      do 215 i=ilower,iupper
        n=(j+1-jlower)*im + i+1-(ilower-1)
        f(n)=f(n)/ac(n)
        ae(n)=ae(n)/ac(n)
        aw(n)=aw(n)/ac(n)
        as(n)=as(n)/ac(n)
        an(n)=an(n)/ac(n)
        ac(n)=1.
  215 continue
  210 continue

      endif

c
c

      itry=0
      
 1122 call umgs2(
     + ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
     + ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
     + ipc,irc,irurd)


      if ((rmax.gt.tol).and.(itry.le.5)) then
        itry=itry+1
        print*,"Retry #",itry
        goto 1122
      endif

      if (rmax.gt.tol) then
        ier = -1
        print*,"Did not converge"
        print*," maximum residual    = ",rmax
        print*," tolerance           = ",tol
      endif

c
c
c
      do 510 j=jlower,jupper
      do 500 i=ilower,iupper
      n=(j+1-jlower)*im + i+1-(ilower-1)
      u(i,j)=q(n)
  500 continue
  510 continue

      return
      end
