      program jacobi
      implicit real*8 (a-h,o-z)
c     /***************************************************************/
c     /* jacobi iteration for poisson problem in 2d, discretized by  */
c     /* second-order central finite-difference on cartesian mesh.   */
c     /* single processor version.                                   */
c     /* m.u. 26/10/2004                                             */
c     /***************************************************************/
c     /* (global) problem size, number of processors, type of solution */
      parameter(m=10,n=10,nproc=1,isol=2)
c
c     /* max. number of iterations allowed & tolerance*/
      parameter(nitmx=100000,tol=1.e-5)
c
c     /* arrays */
      dimension aloc(0:m+1,0:n+1),adebug(0:m+1,0:n+1)
c
c     /* other stuff */
      logical lconverged
      external lconverged
c-----|---------------------------------------------------------------|
c     /* initialize the physical boundaries & set interior to zero*/
      call assign(aloc,m,n,isol)
      call write_array(aloc,m,n,isol,55)
c
      do 100 iter=1,nitmx       !this is the main iteration loop
c     
         call jacobi_step(aloc,adebug,m,n,diffl,valmxl) 
c     
         if(lconverged(diffl,valmxl,tol))goto 110
c     
 100  continue
      write(*,*)'jacobi iteration did not converge! ',iter,
     $     diffl/valmx,tol
 110  continue
      write(*,*)'jacobi iteration converged in ',iter,' steps.'
c
c     /* compute the error of the solution */
      call check_array(aloc,m,n,isol)
c
c     /* write result to file (along with solution) */
      call write_array(aloc,m,n,isol,44)
c
      stop
      end
c----+----------------------------------------------------------------|
      subroutine jacobi_step(a,tmp,m,n,diff,valmx)
      implicit real*8 (a-h,o-z)
c
c     /***************************************************************/
c     /* performs one step of jacobi iteration for 2d matrix in 1d   */
c     /* cartesian domain decomposition along 2nd index.             */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(0:m+1,0:n+1)          matrix in local memory              */
c     /* tmp(0:m+1,0:n+1)        work array                          */
c     /* my_rank                 processor id                        */
c     /*                                                             */
c     /* output:                                                     */
c     /* a(0:m+1,0:n+1)          updated matrix                      */
c     /* diff                    the r.m.s. of the increment         */
c     /*                         relative to the maximum value of the*/
c     /*                         function (for convergence measure)  */
c     /***************************************************************/
c
      dimension a(0:m+1,0:n+1),tmp(0:m+1,0:n+1)
c----+----------------------------------------------------------------|
      pi=4.*atan(1.)
c
      if=1
      il=m
      jf=1
      jl=n
      diff=0.0
      valmx=0.0
      do j=jf,jl
         do i=if,il
            tmp(i,j)=(a(i,j+1)+a(i+1,j)+a(i-1,j)+a(i,j-1))*.25d0
            diff=max(diff,abs(tmp(i,j)-a(i,j)))
            valmx=max(valmx,abs(tmp(i,j)))
         enddo
      enddo
      do j=jf,jl
         do i=if,il
            a(i,j)=tmp(i,j)
         enddo
      enddo
c     
      return
      end subroutine jacobi_step
c-----+----------------------------------------------------------------|
      logical function lconverged(diffl,valmxl,tol)
      implicit real*8 (a-h,o-z)
c     
      lconverged=.false.
c     
      diffl=diffl/valmxl
      if(diffl.le.tol)then
         lconverged=.true.
         write(*,*)' we have converged!',diffl,tol
      endif
c
      return
      end function lconverged
c----+----------------------------------------------------------------|
      subroutine assign(a,m,n,itype)
      implicit real*8 (a-h,o-z)
c
c     /***************************************************************/
c     /* initializes the matrix a to some values.                    */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(m,n)                  matrix in local memory              */
c     /* itype                   type of the requested values.       */
c     /*                         =1: a(i,j)=i*j                      */
c     /*                         =2: a(i,j)=sin(pi*xi)*exp(-pi*yj)   */
c     /*                                                             */
c     /* output:                                                     */
c     /* a(m,n)                  newly assigned matrix               */
c     /***************************************************************/
      dimension a(0:m+1,0:n+1)
c----+----------------------------------------------------------------|
      pi=4.*atan(1.)
c
      if=1
      il=m
      jf=1
      jl=n
      if(itype.eq.2)then
c     /* initialize the interior of the matrix to zero */
         do j=jf,jl
            do i=if,il
               a(i,j)=0.0       !all to zero
            enddo
         enddo
c     /* set physical boundaries:   */
c     /* aij=sin(pi*xi)*exp(-pi*yj) with xi=i/(m+1) \forall 0<=i<=m+1 */
c     /*                            with yj=j/(n+1) \forall 0<=j<=n+1 */
         do j=jf-1,jl+1
            yy=float(j)/float(n+1)
            i=if-1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
            i=il+1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
         enddo
         j=jf-1
         yy=float(j)/float(n+1)
         do i=if-1,il+1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
         enddo
         j=jl+1
         yy=float(j)/float(n+1)
         do i=if-1,il+1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
         enddo
      else
         write(*,*)'assign: this type not implemented: ',itype
         stop
      endif
c     
      return
      end subroutine assign
c----+----------------------------------------------------------------|
      subroutine check_array(a,m,n,itype)
      implicit real*8 (a-h,o-z)
c
c     /***************************************************************/
c     /* computes the error w.r.t. to analytic solution              */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(m,n)                  matrix in local memory              */
c     /* itype                   type of the requested values.       */
c     /*                         =1: a(i,j)=i*j                      */
c     /*                         =2: a(i,j)=sin(pi*xi)*exp(-pi*yj)   */
c     /*                                                             */
c     /* output:                                                     */
c     /***************************************************************/
      dimension a(0:m+1,0:n+1)
c----+----------------------------------------------------------------|
      pi=4.*atan(1.)
c
      if=1
      il=m
      jf=1
      jl=n
      errmx=0.d0
      valmx=0.d0
      if(itype.eq.2)then
         do j=jf,jl
            yy=float(j)/float(n+1)
            do i=if,il
               xx=float(i)/float(m+1)
               sol=solution(xx,yy,itype)
               errmx=max(errmx,abs(a(i,j)-sol))
               valmx=max(valmx,sol)
            enddo
         enddo
         write(*,*)'maximum relative error: ',errmx/valmx,errmx,valmx
      endif
c     
      return
      end subroutine check_array
c----+----------------------------------------------------------------|
      subroutine write_array(a,m,n,itype,iunit)
      implicit real*8 (a-h,o-z)
c
c     /***************************************************************/
c     /* writes 2d array to file, in ascii format with line skip     */
c     /* for gnuplot                                                 */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(m,n)                  matrix in local memory              */
c     /* itype                   type of the requested values.       */
c     /*                         =1: a(i,j)=i*j                      */
c     /*                         =2: a(i,j)=sin(pi*xi)*exp(-pi*yj)   */
c     /* iunit                   unit number for output              */
c     /*                                                             */
c     /* output:                                                     */
c     /***************************************************************/
      dimension a(0:m+1,0:n+1)
c----+----------------------------------------------------------------|
      do j=0,n+1
         yy=float(j)/float(n+1)
         do i=0,m+1
            xx=float(i)/float(m+1)
            sol=solution(xx,yy,itype)
            write(iunit,112)i,j,xx,yy,a(i,j),sol,a(i,j)-sol
         enddo
         write(iunit,111)
      enddo
 111  format(' ')
 112  format(1x,2(i5,1x),10(e15.8,1x))
c
      return
      end subroutine write_array
c----+----------------------------------------------------------------|
      real*8 function solution(x,y,itype)
      implicit real*8 (a-h,o-z)
      pi=4.*atan(1.)
c     
      if(itype.eq.2)then
         solution=sin(pi*x)*exp(-pi*y)
      else
         write(*,*)'solution: this type not implmented! ',itype
         stop
      endif
c     
      return
      end function solution
c----+----------------------------------------------------------------|
