c     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C     %%         Conjugate Gradient Method Solve                 %%
c     %%    -\Delta u=2*pi*pi*sin(pi*x)*sin(pi*y), [0,1]X[0,1]   %%
c     %%    with homogeneous Dirichlet boundary condition        %%
c     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C     r=b-A*x ...... remainder
c     ru,rn ....... means the inner product of (r,r)
c     x ....... CG iterative solution of A*x=b
c     b ....... right hand side matrix
c     p ....... the same defination as the one in CG algorithm
c     ap=A*p .......the same defination as the one in CG algorithm
c     n ...... number of grid points in x-,y-direction
c     h ...... mesh size
c     ex ...... exact solution of the Ax=b
c     it ...... iterative number
c     out.m ...... output file 

      parameter (nn =100 )
      dimension r(0:nn,0:nn),b(0:nn,0:nn),p(0:nn,0:nn),ap(0:nn,0:nn)
      dimension x(0:nn,0:nn)
      pi=2.*asin(1.)
      n=20
      h=1./n
      do j=0,n
        do K=0,n
           xx=j*h
           yy=k*h
           b(j,k)=2.*pi*pi*sin(pi*xx)*sin(pi*yy)*h*h
        enddo
      enddo

C     Set initial guess
      ru=0.
      do j=0,n
        do k=0,n
          x(j,k)=0.
          r(j,k)=b(j,k)
          p(j,k)=b(j,k)
          ru=ru+r(j,k)*r(j,k)
        enddo
      enddo
     
C     Start iteration

      it=0  
  10  it=it+1
      do j=1,n-1
        do k=1,n-1
          ap(j,k)=4.*p(j,k)-p(j-1,k)-p(j+1,k)-p(j,k+1)-p(j,k-1)
        enddo
      enddo
      
      apa=0.
      do j=1,n-1
        do k=1,n-1
          apa=apa+p(j,k)*ap(j,k)
        enddo
      enddo
      
      rn=0.
      alpha=ru/apa
      do j=1,n-1
        do k=1,n-1
          x(j,k)=x(j,k)+alpha*p(j,k)
          r(j,k)=r(j,k)-alpha*ap(j,k)
          rn=rn+r(j,k)*r(j,k)
        enddo
      enddo
      
      beta=rn/ru
      do j=1,n-1
        do k=1,n-1
          p(j,k)=r(j,k)+beta*p(j,k)
        enddo
      enddo  
      ru=rn
      write (*,*) 'iterative number=',it,'remainder=',sqrt(rn)

      if(sqrt(rn).gt.1.e-8 ) go to 10

      open(11,file='out.m',status='unknown')
      write(11,*) 'N=',n,';'
      write(11,*) 'sol=['
      do j=0,n
        do k=0,n
          xx=j*h
          yy=k*h
          ex=sin(pi*xx)*sin(pi*yy)
          write(11,*) xx,yy,x(j,k),ex,abs(x(j,k)-ex)
        enddo
      enddo
      write(11,*) '];' 
      close(11)
      
      stop
      end