c     Main
      implicit real*8 (a-h, o-z)
      parameter(mesh=201, nx=mesh, ny=mesh)
      common/ivect/ivect, /itest/  itest
      ivect=0
      itest=0
      call initia
      call micgs(1d-14,n)
      call second
      call clock(ctime,3)
      call micgs(1d-14,n)
      call clock(ctime,5)
      write(*,'('' ivect='',i1,3x,''LDU='',i1)')ivect,1-itest
      write(*,'('' N='',i4,3x,''Loop='',i5,3x,''CPU time='',1pe10.3)')
     *      mesh,n,ctime 
      stop
      end
c
c     patch for clock routine when running on Unix
c
      subroutine clock(ctime, n)
      implicit real*8 (a-h, o-z)
      real*4 timear(2), dtime
      if (n.eq.3) then
        ctime=dtime(timear)
        ctime=0d0
      end if
      if (n.eq.5) ctime=dtime(timear)+ctime
      return
      entry xclock(remain, n)
      remain=1e10
      return
      end  
c
c *** initial setting ***
c
      subroutine initia
      implicit real*8 (a-h, o-z)
      parameter(mesh=201, nx=mesh, ny=mesh)
      parameter(im=nx, jm=ny, im1=im-1, jm1=jm-1, pi=3.14159265d0)
      parameter(iz=nx-2, ir=ny-2)
      common /scalar/ rho(nx,ny), psi(nx,ny)
      common /poissn/ ca(iz),cb(iz),cc(iz),cd(ir),ce(ir),cf(ir)
      dimension r(10000), dr(10000)
c
      dx=1d0/float(ir+1)
      do 100 j=1, ir+1
      r(j)=dx*(float(j)-0.5d0)
  100 dr(j)=dx
c
c     coefficients of poisson equation
c      (d/dz)**2
      do 300 i=2, iz
        ca(i)=1d0/dx**2
        cb(i)=2d0/dx**2
        cc(i)=1d0/dx**2
  300 continue
      ca(1)=1d0/dx**2
      cb(1)=2d0/dx**2
      cc(1)=ca(1)
c
      do 310 j=2, ir
        cd(j)=(2*r(j)+dr(j))/
     *           ((r(j)+dr(j)/4-dr(j-1)/4)*dr(j)*(dr(j-1)+dr(j)))
        ce(j)=2*r(j)/((r(j)+dr(j)/4-dr(j-1)/4)*dr(j)*dr(j-1))
        cf(j)=(2*r(j)-dr(j-1))/
     *           ((r(j)+dr(j)/4-dr(j-1)/4)*dr(j-1)*(dr(j-1)+dr(j)))
  310 continue
c
      cd(1)=4.0/(dr(1)*(2*r(1)+dr(1)))
      ce(1)=cd(1)
c
c
c  ** setting **
c     for test
      do 200 j=1, ir+1
      do 200 i=1, iz+1
        rho(i,j)=1d0
  200 continue
      return
      entry second
      do 210 j=1, ir+1
        if(r(j).le.0.5d0) then
          do 220 i=1, nx
            rho(i,j)=1d0
  220     continue
        else
          do 230 i=1, nx
            rho(i,j)=0d0
  230     continue
        end if
  210 continue
      return
      end
c
c     preconditioning incomplete lu decomposition
c     conjugate gradient squared method
c
      subroutine micgs  (eps2, n)
      implicit real*8(a-h,o-z)
      parameter (mesh=201, nx=mesh, ny=mesh)
      parameter (iz=nx-2, ir=ny-2, iz1=iz+1, ir1=ir+1, izr=iz*ir)
      common /scalar/ rho(nx,ny), psi(nx,ny)
      common /m1/ x(1-iz:izr+iz), y(1-iz:izr+iz), w(1-iz:izr+iz)
      common /m2/ r0(1-iz:izr+iz), r(1-iz:izr+iz), b(1-iz:izr+iz),
     *          p(1-iz:izr+iz), e(izr), h(izr)
      common /A/a(1:izr), ub(1-iz:izr), uc(1-iz:izr),
     *          db(1:izr), dc(1:izr), d(1-iz:izr),
     *          ud(1-iz:izr), dd(0:izr)
      common /poissn/ ca(iz), cb(iz), cc(iz), cd(ir), ce(ir), cf(ir)
      common /itest/  itest
      data initia/1/
c      itest=1
      if(initia.ne.1) goto 1000
      initia=0
c
      do 10 j=1, ir
      do 10 i=1, iz
        k=iz*(j-1)+i
        a(k) =  cb(i)+ce(j)
        ub(k) = -ca(i)
        db(k) = -cc(i)
        uc(k) = -cd(j)
        dc(k) = -cf(j)
   10 continue
c
      do 12 k=1-iz, -1
        ud(k)=0d0
   12 continue
      do 14 k=0, izr
        ud(k)=0d0
        dd(k)=0d0
   14 continue
c
      do 20 j=1, ir
c
c       upper boundary
        k=iz*j
        ub(k)=0d0
        dd(k)=-ca(iz)
c
c       lower boundary
        k=iz*j-iz+1
        db(k)=0d0
        ud(k)=-cc(1)
   20 continue
c
      do 30 i=1, iz
c
c       outer boundary
        k=iz*(ir-1)+i
        uc(k)=0d0
c
c       inner boundary
        dc(i)=0d0
   30 continue
c
      do 40 i=1-iz, izr+iz
        x(i)=0d0
        y(i)=0d0
        w(i)=0d0
        b(i)=0d0
        p(i)=0d0
   40 continue
c
 1000 continue
c
c     source term
c
      call boundp
c
      do 100 j=1, ir
        do 90 i=1, iz
          k=iz*(j-1)+i
          z=rho(i+1,j+1)
          b(k)=-(z+sign(z,z-1d0))/2
          x(k)=psi(i+1,j+1)
   90   continue
c
c       upper boundary
*       b(iz*j)=b(iz*j)+ca(iz)*psi(iz1+1,j+1)
  100 continue
c
c       outer boundary
      do 110 i=1, iz
        b(iz*(ir-1)+i)=b(iz*(ir-1)+i)+cd(ir)*psi(i+1,ir1+1)
  110 continue
c
c     b = (ldu)**(-1) b
c
      if(itest.ne.1) call mulldu(b)
c
c     y = a x(0)
c
      call multia(x, y)
c
c     r(0) = b - a x(0) ; p(0) = e(0) = r(0)
c
      do 200 i=1, izr
        r0(i)=b(i)-y(i)
        p(i)=r0(i)
        e(i)=r0(i)
        r(i)=r0(i)
  200 continue
      n=0
      b2=prd(b,b)
c
 2000 continue
c
c     y = a p
c
      call multia(p, y)
c
c     alpha(n) = (r(0), r(n))/(r(0), a p(n))
c
      r0rn=prd(r0,r)
      alpha=r0rn/prd(r0,y)
c
c     h(n+1) = e(n) - alpha(n) a p(n)
c     w = e(n) + h(n+1)
c
      do 300 i=1, izr
        h(i)=(e(i)-alpha*y(i))
        w(i)=(e(i)-alpha*y(i))+e(i)
  300 continue
c
c     y = a ( e(n) + h(n+1) )
c
      call multia(w, y)
c
c     r(n+1) = r(n) - alpha(n) a ( e(n) + h(n+1) )
c     x(n+1) = x(n) + alpha(n) ( e(n) + h(n+1) )
c
      do 400 i=1, izr
        r(i)=r(i)-alpha*y(i)
        x(i)=x(i)+alpha*w(i)
  400 continue
c
c     beta(n) = (r(0),r(n+1))/(r(0),r(n))
c
      beta=prd(r0,r)/r0rn
c
c     e(n+1) = r(n+1) + beta(n) h(n+1)
c     p(n+1) = e(n+1) + beta(n) (h(n+1) + beta(n) p(n))
c
      do 500 i=1, izr
        e(i)=(r(i)+beta*h(i))
        p(i)=(r(i)+beta*h(i))+beta*(h(i)+beta*p(i))
  500 continue
c
      n=n+1
c
c     until |r(n)| < eps |b| do.
c
      r2=prd(r,r)
      if(r2.gt.eps2*b2) goto 2000
c
c     return to psi
c
      do 700 j=1, ir
      do 700 i=1, iz
        psi(i+1,j+1)=x(iz*(j-1)+i)
  700 continue
c
      return
      end
c
c     y = a x
c
      subroutine multia (x, y)
      implicit real*8(a-h,o-z)
      parameter (mesh=201, nx=mesh, ny=mesh)
      parameter (iz=nx-2, ir=ny-2, iz1=iz+1, ir1=ir+1, izr=iz*ir)
      common /A/a(1:izr), ub(1-iz:izr), uc(1-iz:izr),
     *          db(1:izr), dc(1:izr), d(1-iz:izr),
     *          ud(1-iz:izr), dd(0:izr)
      dimension x(-iz+1:izr+iz), y(-iz+1:izr+iz)
      common /itest/  itest
c
      do 100 i=1, izr
        y(i)=a(i)*x(i)+ub(i)*x(i+1)+db(i)*x(i-1)
     *         +uc(i)*x(i+iz)+dc(i)*x(i-iz)
     *         +ud(i)*x(i+iz-1)+dd(i)*x(i-iz+1)
  100 continue
c
c     itest=1 means no lu
c
      if (itest.eq.1) return
c
c     y = (ldu)**(-1) y
c
      call mulldu(y)
c
      return
      end
c
c     y = (ldu)**(-1) y
c
      subroutine mulldu (y)
      implicit real*8(a-h,o-z)
      parameter (mesh=201, nx=mesh, ny=mesh)
      parameter (iz=nx-2, ir=ny-2, iz1=iz+1, ir1=ir+1, izr=iz*ir)
      common/A/ a(1:izr), ub(1-iz:izr), uc(1-iz:izr),
     *          db(1:izr), dc(1:izr), d(1-iz:izr),
     *          ud(1-iz:izr), dd(0:izr)
      common/ivect/ivect
      dimension y(-iz+1:izr+iz)
      dimension ln(0:iz+ir-1), l(1:izr)
      dimension parm(7)
      data   parm/0.95d0, 0.90d0, 0.75d0, 0.5d0, 0d0, -1d0, -3d0/
      data   initia/1/
c
c     initialize
c
      if (initia.ne.1) goto 1000
c
c     incomplete ldu
c
      do 1 i=1-iz, 0
        d(i)=0d0
    1 continue
      do 3 k=1, 7
        u=parm(k)
        nparm=0
c
        do 2 i=1, izr
          w=a(i)-db(i)*d(i-1)*ub(i-1)-dc(i)*d(i-iz)*uc(i-iz)
     *        -dd(i)*d(i-iz+1)*ud(i-iz+1)
     *        -u*(db(i)*d(i-1)*uc(i-1)+dc(i)*d(i-iz)*ub(i-iz)
     *           +dd(i)*d(i-iz+1)*ub(i-iz+1)+db(i)*d(i-1)*ud(i-1))
          if(w.lt.a(i)*1e-1) nparm=1
          d(i)=1d0/w
    2   continue
        if(nparm.eq.0) goto 5
    3 continue
    5 continue
c
c     prepare for zensin, koutai dainyu
c
      if(ivect.eq.1) then
        ln(0)=0
        k=0
        do 20 m=1, iz
        do 10 i=m, 1, -1
          j=m+1-i
          if(j.le.ir) then
            k=k+1
            l(k)=iz*(j-1)+i
          end if
   10   continue
        ln(m)=k
   20 continue
c
      do 40 m=iz+1, iz+ir-1
        do 30 i=iz, m-ir+1, -1
          j=m+1-i
          k=k+1
          l(k)=iz*(j-1)+i
   30   continue
        ln(m)=k
   40 continue
      end if
c
      initia=0
c
 1000 continue
c
      if(ivect.eq.1) then
c
c     << vector >>
c     y = l**(-1) y
c
      do 100 j=1, iz+ir-1
*voption indep(y)
*vocl loop,novrec(y)
      do 100 i=ln(j-1)+1, ln(j)
        y(l(i)) = d(l(i)) * ( y(l(i)) - dc(l(i))*y(l(i)-iz)
     *                        - db(l(i))*y(l(i)-1)
     *                        - dd(l(i))*y(l(i)-iz+1) )
  100 continue
c
c     y = (du)**(-1) y
c
      do 200 j=iz+ir-1, 1, -1
*voption indep(y)
*vocl loop,novrec(y)
      do 200 i=ln(j-1)+1, ln(j)
        y(l(i)) = y(l(i)) - d(l(i)) * ( ub(l(i))*y(l(i)+1)
     *                                  + uc(l(i))*y(l(i)+iz)
     *                                  + ud(l(i))*y(l(i)+iz-1) )
  200 continue
c
      else
c
c     << scalar >>
c     y = l**(-1) y
c
      do 1100 i=1, izr
        y(i) = d(i) * ( y(i) - dc(i)*y(i-iz) - db(i)*y(i-1)
     *                       - dd(i)*y(i-iz+1) )
 1100 continue
c
c     y = (du)**(-1) y
c
      do 1200 i=izr-1, 1, -1
        y(i) = y(i) - d(i) * ( ub(i)*y(i+1) + uc(i)*y(i+iz)
     *                       + ud(i)*y(i+iz-1) )
 1200 continue
c
      end if
c
      return
      end
c
c     inner product
c
      real*8 function prd(a,b)
      implicit real*8(a-h,o-z)
      parameter (mesh=201, nx=mesh, ny=mesh)
      parameter (iz=nx-2, ir=ny-2, iz1=iz+1, ir1=ir+1, izr=iz*ir)
      dimension a(1-iz:izr+iz), b(1-iz:izr+iz)
c
      prd=0d0
c
      do 100 i=1, izr
        prd=prd+a(i)*b(i)
  100 continue
c
      return
      end
c
c   program for boundary condition to gravitational potential
c
      subroutine boundp
      implicit real*8(a-h,o-z)
      parameter (mesh=201, nx=mesh, ny=mesh, im=nx, jm=ny)
      parameter (iz=nx-2, ir=ny-2, iz1=iz+1, ir1=ir+1)
      common /scalar/ rho(nx,ny),  psi(nx,ny)
      common /mesh  / xb(nx), yb(ny),
     *                xf(nx), yf(ny),
     *                dxb(nx), dyb(ny),
     *                dxf(nx), dyf(ny), dvol(ny)
c
c  boundary at i=ir (right)
c
      do 100 i=2, im
        psi(i,jm)=1d0
  100 continue
c
c  boundary at i=iz (upper)
c
*     do 200 j=2, jm
*       psi(im,j)=0.0
* 200 continue
c
      return
      end