Line Jacobi (Fortran Version)

      program ljacobi
      parameter(mmax=50,nmax=50,mx=mmax+1,ny=nmax+1)
      parameter(tol=.005)
      dimension u(0:mx,0:ny)
      dimension v(0:mx,0:ny)
      dimension e(0:mx,0:ny)
      dimension a(nmax), b(nmax), c(nmax), d(nmax)
c*********************************************
c
c    Line Jacobi algorithm for
c    Laplace difference equation with
c    Dirichlet boundary conditions
c
c**********************************************
c    Initialize  arrays and set boundary conditions
c**********************************************
      dx = 1./mx
      do 2 m = 1,mmax
      x = m*dx
      v(m,0) = x*x 
      v(m,ny) = x*x - 1.
2     continue
      dy = 1./ny
      do 3 n = 1,nmax
      y = n*dy
      v(0,n) = -y*y
      v(mx,n) = 1. - y*y
3     continue
      do 99 j = 1,nmax
      do 99 i = 1,mmax
      x = i*dx
      y = j*dy
      e(i,j) = x*x - y*y
      v(i,j) = 0.0
99    continue
 
c********************************************
c   start Line Jacobi iterations
c********************************************
      do 77 k = 1,100
 
c**********************************************
c     Perform 10 Line Jacobi iterations and check error
c**********************************************
c  step across the m columns      
      do 6 m  = 1,mmax
c  build tridiagonal system for column m
      do 55 n = 1,nmax
      a(n) = -1.
      b(n) = 4.
      c(n) = -1.
      d(n) = v(m-1,n) + v(m+1,n)
55    continue
      d(1) = d(1) + v(m,0)
      d(nmax) = d(nmax) + v(m,ny)
c   start forward substitution
      do 10 i = 2,nmax
      ratio = a(i)/b(i-1)
      b(i) = b(i) - ratio*c(i-1)
      d(i) = d(i) - ratio*d(i-1)
10    continue
c   initialize backward substitution
      d(nmax) = d(nmax)/b(nmax)
      u(m,nmax) = d(nmax)
c   finish bacward substitution
      do 20 i = nmax-1,1,-1
      d(i) = (d(i) - c(i)*d(i+1))/b(i)
      u(m,i) = d(i)
20    continue
6     continue
 
c   write u back into v to prepare for next iteration
      do 7 n = 1,nmax
      do 7 m = 1,mmax
      v(m,n) = u(m,n) 
7     continue
 
 
c   output residual and error information
      residual = 0.0
      emax = 0.0
      do 88 n = 1,mmax
      do 88 m = 1,mmax
      res = 4.*v(m,n) - (v(m-1,n)+v(m,n-1)+v(m+1,n)+v(m,n+1))
      residual = amax1(abs(res),residual)
      error = v(m,n) - e(m,n)
      emax = amax1(abs(error),emax)
88    continue
      write(6,200) k,residual,emax
      if(emax.lt.tol) goto 999
 
77    continue
999   write(6,300)
 
200   format('at iteration ',i5' residual = ',e10.2' error = ',e10.2)
 
300   format('done')
      end