C********************************************************************
C--- HOMEWORK #6
C--- Programmed by Sung-Hoon Ko(shko@npac.syr.edu)
C---=================================================================
C--- Numericl Solution to Partial Differential Equations.
C-
C--- Solve Laplace's equation in 2 dimension using the 
C--- Conjugate Gradients and Jacobi, and display the error
C--------------------------------------------------------------------
C-
C-  We have a Laplace's matrix with boundary values 1 and
C-  interior values 0.
C-  We use Conjugate Gradients and Jacobi to solve Laplace's equation.
C-  However, Conjugate Gradients use matrix A and b 
C-  to solve equation Ax=b.
C-  How can we get matrix A and b from Laplace's matrix?
C-  Suppose we have 5*5 Laplace's matrix.
C-  
C-       b01  b02  b03  b04  b05
C-       b06  x01  x02  x03  b07
C-       b08  x04  x05  x06  b09
C-       b10  x07  x08  x09  b11
C-       b12  b13  b14  b15  b16
C-  
C-     where b is boundary values.
C-  Then   x01 = (x02+x04+b02+b06)/4
C-         x02 = (x01+x03+x05+b03)/4
C-  and so on.
C-
C-  We can treat it as a system of 9 equations in 9 unknowns
C-  with a 9-by-9 coefficient matrix:
C-
C-  +--                         --+   +-   -+     +-       -+
C-  |  4 -1  0 -1  0  0  0  0  0  |   | x01 |     | b02+b06 |
C-  | -1  4 -1  0 -1  0  0  0  0  |   | x02 |     | b03     |
C-  |  0 -1  4  0  0 -1  0  0  0  |   | x03 |     | b03+b07 |
C-  | -1  0  0  4 -1  0 -1  0  0  |   | x04 |     | b08     |
C-  |  0 -1  0 -1  4 -1  0 -1  0  |   | x05 |  =  |   0     |
C-  |  0  0 -1  0 -1  4  0  0 -1  |   | x06 |     | b09     |
C-  |  0  0  0 -1  0  0  4 -1  0  |   | x07 |     | b10+b13 |
C-  |  0  0  0  0 -1  0 -1  4 -1  |   | x08 |     | b14     |
C-  |  0  0  0  0  0 -1  0 -1  4  |   | x09 |     | b11+b15 |
C-  +--                         --+   +-   -+     +-       -+
C-               A                       x     =       b
C-
C-  Subroutine get_A and get_b find this matrix from 
C-  Laplace's matrix then use subroutine ConGrad to
C-  perform Conjugate Gradients.
C-
C*********************************************************************

C---  Note: Fortran90 version

      PROGRAM  HW6_LAPLACE_ConGrad_Jacobi
      IMPLICIT NONE
      INTEGER, PARAMETER :: M=27, N=M-2
      REAL    Cx(1:N*N,1), Jx(M,M)
      REAL    Ceps(100), Jeps(100)   ! store the errors in array to draw
      INTEGER  Citer, Jiter
      EXTERNAL dp_initialize  ! External routines to initialize
      EXTERNAL dp_draw_graph  ! and draw

      CALL  Start_ConGrad(Cx, Ceps, Citer)

      CALL  Start_Jacobi(Jx, Jeps, Jiter)

      WRITE(6,11) ' *** Number of iterations ***'
      WRITE(6,33)'Conjugate Gradients :', Citer
      WRITE(6,33)'Jacobi :', Jiter
 11   FORMAT(///,A,/)
 33   FORMAT(8X,A,I3,/)
      
      CALL  dp_initialize
      CALL  dp_draw_graph(Ceps, Jeps, Citer, Jiter)
      PAUSE
           
      STOP
      END



C********************************************************************
C---
      SUBROUTINE Start_ConGrad(x, eps, iter)
      IMPLICIT NONE
      INTEGER, PARAMETER :: M=27, N=M-2
      REAL    A(1:N*N, 1:N*N)
      REAL    x(1:N*N,1), b(1:N*N,1)
      REAL    eps(*)   ! store the errors in array to draw
      INTEGER iter

C--- Find coefficient matrix A
      CALL  get_A(A)

C--- Find matrix b
      CALL  get_b(b)

C--- Perform Conjugate Gradients
      CALL  ConGrad(A,x,b,eps,iter)

      RETURN
      END



C=====================================================================
C---
C---  Perform Conjugate Gradients on equation Ax=b 
C---  where A is symmetric positive definite.
C---  Returns updated functions, the number of total iterations,
C---  writes error estimate in maxError

      SUBROUTINE ConGrad(A,x,b, Error, k)
      IMPLICIT NONE
      INTEGER, PARAMETER :: M=27, N=M-2
      REAL  A(1:N*N, 1:N*N)
      REAL, ARRAY(1:N*N,1) :: x, b
      REAL  Error(50), maxError
      INTEGER k
      REAL, ARRAY(1:N*N,1) :: rnew, rold, p, xnew
      REAL  beta, alpha
      REAL  tbeta(1,1), talpha(1,1)  ! temporary variables

      x=0.
      rnew=b
      k=0
      maxError=1.0
      DO WHILE(maxError .GT. 0.01)
         k=k+1
         IF(k .EQ. 1)  THEN
             p = rnew
         ELSE
             tbeta=matmul(transpose(rnew),rnew)/
     *             matmul(transpose(rold),rold)
             beta=tbeta(1,1)
             p=rnew+beta*p

         ENDIF         
         talpha=matmul(transpose(rnew),rnew)/
     *          matmul(matmul(transpose(p),A),p)
         alpha=talpha(1,1)
         xnew = x + alpha*p
         rold = rnew
         rnew = rnew - alpha*matmul(A,p)
         maxError=maxval(abs(xnew-x))
         Error(k)=maxError
         x = xnew
      ENDDO

      RETURN
      END



C------------------------------------------------------
C--- Find coefficient matrix A from Laplace's matrix
      
      SUBROUTINE  get_A(A)
      INTEGER, PARAMETER :: N=25
      REAL    A(1:N*N, 1:N*N)
      INTEGER Atemp(1:N, 1:N)
      INTEGER i,j,k

      A = 0
C--- Put diagonal 4
      FORALL(i=1:N*N, j=1:N*N, (i .EQ. j))
     *    A(i,j) = 4.
      Atemp = RESHAPE( [1:N*N], [N,N] )

      DO i=1,N
        DO j=1,N
          k = Atemp(j,i)
          IF(j .GT. 1)  A(k,k-1)=-1.
          IF(j .LT. N)  A(k,k+1)=-1.
          IF(i .GT. 1)  A(k,k-N)=-1.
          IF(i .LT. N)  A(k,k+N)=-1.
        ENDDO
      ENDDO
      
      RETURN
      END


C-------------------------------------------------
C--- Find matrix b from Laplace's matrix

      SUBROUTINE  get_b(b)
      INTEGER, PARAMETER :: M=27, N=M-2
      REAL    b(1:N*N,1)
      REAL    btemp(1:M,1:M)
      INTEGER i,j,k

C---  Install boundary values
      btemp(1,1:M)=1
      btemp(1:M,1)=1
      btemp(M,1:M)=1
      btemp(1:M,M)=1

C--- Initialize interior points
      btemp(2:M-1,2:M-1)=0

C--- Then add four directions, which is same as
C--- value of matrix b
      btemp(2:M-1,2:M-1) = btemp(1:M-2,2:M-1) 
     *                   + btemp(3:M,2:M-1)
     *                   + btemp(2:M-1,1:M-2) 
     *                   + btemp(2:M-1,3:M) 
      
C--- Put the values of b
      DO i=2,M-1
        DO j=2,M-1
          k=(i-2)*(M-2)+j-1
          b(k,1) =  btemp(i,j)
        ENDDO
      ENDDO

      RETURN 
      END

 
C********************************************************************
C--- HOMEWORK #5.(It was homework 5)
C--- Numericl Solution to Partial Differential Equations.
C---
C--- Solve Laplace's equation in 2 dimension using the 
C--- point-Jacobi iteration.
C--------------------------------------------------------------------
C--- Programmed by Sung-Hoon Ko(shko@npac.syr.edu)
C--- (It was modified for homework 6)
C********************************************************************


      SUBROUTINE  Start_JACOBI(GRID, eps, iter)
      IMPLICIT NONE
      REAL    GRID(1:27,1:27), eps(*)
      REAL    error
      INTEGER iter,i,j

C---  Install boundary values
      GRID(1,1:27)=1
      GRID(1:27,1)=1
      GRID(27,1:27)=1
      GRID(1:27,27)=1

C--- Initialize interior points
      GRID(2:26,2:26)=0

      iter=0
      error=1.0
      DO WHILE(error .GT. 0.01)
         CALL JACOBI(GRID, error)
         iter=iter+1
         eps(iter)=error
      ENDDO

      RETURN
      END



C=====================================================================
C---
C---  Perform Point-Jocobi on grid points G.
C---  Returns updated functions, writes error estimate in maxError.

      SUBROUTINE JACOBI(G, MaxError)
      REAL    G(1:27,1:27)   , Gnew(1:27,1:27)
      REAL    MaxError

      Gnew = G
      Gnew(2:26,2:26) = 0.25*(  G(1:25,2:26) 
     *                        + G(3:27,2:26)
     *                        + G(2:26,1:25) 
     *                        + G(2:26,3:27) )
      
C----Find the maximum error 
      maxError=maxval(abs(Gnew-G))     

      G = Gnew    ! Install updated Grid values

      RETURN
      END
      




















