c  
c  solution of Laplace's eq on rectangle
c           0 < x < xmax, 0 < y < ymax
c  
c  boundary conditions specified in function
c    statements at end of code
c  
c  CGM used to solve linear system
c  
      implicit real*8(a-h,o-z)
      parameter(nx=30,ny=30,n=900,tol=0.00001)
      parameter(xmax=1.0,ymax=1.0)
      dimension b(n),v(n),w(n),q(n),r(n),d(n)
      dx=xmax/(nx+1)
      dy=ymax/(ny+1)
      bb=dx*dx/(dy*dy)
c  
c  define starting value for v
c  
      do 4 i=1,n
        v(i)=0.0
    4   continue
      call right(n,nx,ny,dx,dy,bb,b)
      call mult(n,nx,ny,bb,v,q)
      rr=0.0
      do 8 i=1,n
        w(i)=v(i)
        r(i)=b(i)-q(i)
	d(i)=r(i)
    8   rr=rr+r(i)**2
c  
c  CGM iteration
c  
      do 20 it=1,500
        if(it.eq.1) beta=0.0
        if(it.ne.1) beta=rr/rr0
	do 10 i=1,n
   10	  d(i)=r(i)+beta*d(i)
        call mult(n,nx,ny,bb,d,q)
	dAd=0.0
	do 12 i=1,n
   12	  dAd=dAd+d(i)*q(i)
        alpha=rr/dAd
	rr0=rr
	rr=0.0
        do 14 i=1,n
	  v(i)=v(i)+alpha*d(i)
	  r(i)=r(i)-alpha*q(i)
   14	  rr=rr+r(i)**2
	err1=0.0
	do 16 i=1,n
	  err1=err1+(v(i)-w(i))**2
	  w(i)=v(i)
   16   continue
        if(sqrt(err1).lt.tol) go to 22
   20   continue
c  
c  format statements
c  
   22 open(unit=6,file="data")
      write(6,103)
  103 format('U:=[')
      x0=0.0
      y0=0.0
      write(6,104)
      write(6,100) (ix*dx,y0,gB(ix*dx),ix=0,nx)
      write(6,107) xmax,y0,gB(xmax)
      write(6,105)
      do 30 iy=1,ny
        write(6,104)
        write(6,100) x0,iy*dy,gL(iy*dy)
        write(6,100) (ix*dx,iy*dy,v((iy-1)*nx+ix),ix=1,nx)
        write(6,107) xmax,iy*dy,gR(iy*dy)
   30   write(6,105)
      write(6,104)
      write(6,100) (ix*dx,ymax,gT(ix*dx),ix=0,nx)
      write(6,107) xmax,ymax,gT(xmax)
      write(6,106)
  100 format('[',f8.4,',',f8.4,',',e12.4,'],')
  107 format('[',f8.4,',',f8.4,',',e12.4,']')
  104 format('[')
  105 format('],')
  106 format(']];')
      close(unit=6)
      open(unit=7,file="iterations")
      write(7,201)  it
  201 format(i6)
      close(unit=7)
      end
  
c  
c  calculate q = Ap
c  
      subroutine mult(n,nx,ny,bb,p,q)
      implicit real*8(a-h,o-z)
      dimension q(n),p(n)
      np=n-nx
      do 10 i=1,n
        q(i)=2.0*(1.0+bb)*p(i)
        if(i.gt.1) q(i)=q(i)-p(i-1)
        if(i.lt.n) q(i)=q(i)-p(i+1)
        if(i.le.np) q(i)=q(i)-bb*p(i+nx)
        if(i.gt.nx) q(i)=q(i)-bb*p(i-nx)
   10   continue
      do 20 j=1,ny-1
        i=j*nx
        q(i)=q(i)+p(i+1)
   20   q(i+1)=q(i+1)+p(i)	
      return
      end
c  
c  calculate b
c  
      subroutine right(n,nx,ny,dx,dy,bb,b)
      implicit real*8(a-h,o-z)
      dimension b(n)
      do 1 i=1,n
    1   b(i)=0.0 
      do 2 ix=1,nx
        x=ix*dx
        b(ix)=b(ix)+bb*gB(x)
	ixx=nx*(ny-1)+ix
    2	b(ixx)=b(ixx)+bb*gT(x)
      do 4 iy=1,ny
        y=iy*dy
	iyy=nx*(iy-1)
	b(iyy+1)=b(iyy+1)+gL(y)
    4   b(iyy+nx)=b(iyy+nx)+gR(y)
      return
      end     
c  
c  specify boundary conditions
c      
      function gT(x)
      implicit real*8(a-h,o-z)
c      gT=1.0
      gT=sin(4*x*3.14159)
      return
      end
      
      function gB(x)
      implicit real*8(a-h,o-z)
      gB=-2*x*(1-x**4)
      return
      end
      
      function gR(y)
      implicit real*8(a-h,o-z)
      gR=0.0
      return
      end
      
      function gL(y)
      implicit real*8(a-h,o-z)
      gL=0.0
      return
      end