PROGRAM ELLIPTIC_SOLVER
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************            VARIABLES AND PARAMETERS            *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

INTEGER(I4B)::i,j,k,l
INTEGER(I4B),PARAMETER::N=21,MAXITRS=20,maxw=501
REAL(DP)::DX,DY,w,ERROR
REAL(DP), PARAMETER::deltaw=0.002
REAL(DP),DIMENSION(:,:),ALLOCATABLE::X,Y,S,U0,UOLD,UNEW,R
ALLOCATE(X(N,N),Y(N,N),S(N,N),U0(N,N),UOLD(N,N),UNEW(N,N),R(N,N))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************                   GRID SIZING                  *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

CALL GRID(X,Y,DX,DY,N)
    WRITE(*,*)'The Matrix X(i,j) is given as:'
    DO i=LBOUND(X,1),UBOUND(X,1)
        WRITE(*,*) (X(i,j), j=LBOUND(X,2),UBOUND(X,2))
    END DO
    
    WRITE(*,*)'The Matrix Y(i,j) is given as:'
    DO i=LBOUND(Y,1),UBOUND(Y,1)
        WRITE(*,*) (Y(i,j), j=LBOUND(Y,2),UBOUND(Y,2))
    END DO
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************               ANALYTICAL SOLUTION              *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO j=1,N
    DO i=1,N
        U0(i,j)=(X(i,j))**2-(Y(i,j))**2
    END DO
END DO

    WRITE(*,*)'The Analytical Solution of U0(i,j) is given as:'
    DO i=LBOUND(U0,1),UBOUND(U0,1)
        WRITE(*,*)(U0(i,j), j=LBOUND(U0,2),UBOUND(U0,2))
    END DO
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************         INITIAL AND BOUNDARY CONDITIONS        *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO i=1,N
    S(i,1)=U0(i,1)              !Boundary Condition from Analytical Solution
    S(i,N)=U0(i,N)              !Boundary Condition from Analytical Solution
END DO

DO j=1,N
    S(1,j)=U0(1,j)              !Boundary Condition from Analytical Solution
    S(N,j)=U0(N,j)              !Boundary Condition from Analytical Solution
END DO

DO i=2,N-1
    DO j=2,N-1
        S(i,j)=0.0                               !Initial-Value from Problem
    END DO
END DO

DO i=1,N
    DO j=1,N
        UOLD(i,j)=S(i,j)           !Equating Values from S(i,j) to UOLD(i,j)
    END DO
END DO

    WRITE(*,*)'The Initial Value Matrix UOLD(i,j) is given as:'
    DO i=LBOUND(UOLD,1),UBOUND(UOLD,1)
      WRITE(*,*) (UOLD(i,j), j=LBOUND(UOLD,2),UBOUND(UOLD,2))
    END DO
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************               NUMERICAL SOLUTION               *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO i=1,N                                !Boundary Conditions Remain Constant
    UNEW(i,1) = UOLD(i,1)
    UNEW(i,N) = UOLD(i,N)
END DO

DO j=1,N        
    UNEW(1,j) = UOLD(1,j)
    UNEW(N,j) = UOLD(N,j)
END DO

UNEW(i,j)=UOLD(i,j)
ERROR=0.0
w=1.0
DO l=1,maxw
    DO k=1,MAXITRS
        DO j=2,N-1
            DO i=2,N-1
                CALL PSOR(UOLD,UNEW)
                    CALL ERROR_CALCULATION(UOLD,UNEW,ERROR)
                        CALL RR(UNEW,U0,R)
                           
                UOLD(i,j)=UNEW(i,j)
            END DO
        END DO
         WRITE(*,*)'The Initial Value Matrix R(i,j) is given as:'
                            DO i=LBOUND(R,1),UBOUND(R,1)
                            WRITE(*,*) (R(i,j), j=LBOUND(R,2),UBOUND(R,2))
                            END DO
    END DO
    WRITE(*,*)w,ERROR
    ERROR=0.0
    w=w+deltaw
END DO

CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************                   SUBROUTINES                  *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE GRID(X,Y,DX,DY,N)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    INTEGER(I4B)::i,j
    REAL(DP), DIMENSION(N,N)::X,Y
    INTEGER(I4B),INTENT(IN)::N
    REAL(DP),INTENT(OUT)::DX,DY
    
    DX=1.0/(FLOAT(N)-1.0)
    DY=1.0/(FLOAT(N)-1.0)
    
    DO j=1,N
        DO i=1,N
            X(i,j)=(i-1)*DX
            Y(i,j)=(j-1)*DY
        END DO
    END DO
END SUBROUTINE GRID

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE PSOR(UOLD,UNEW)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP)::SUM1
    REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD
    REAL(DP),DIMENSION(N,N),INTENT(OUT)::UNEW

SUM1=0.25*(UOLD(i+1,j)+UNEW(i-1,j)+UOLD(i,j+1)+UNEW(i,j-1)-(4.0*UOLD(i,j)))
UNEW(i,j)=UOLD(i,j)+(w*SUM1)
END SUBROUTINE PSOR

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE ERROR_CALCULATION(UOLD,UNEW,ERROR)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD,UNEW
    REAL(DP),INTENT(OUT)::ERROR

    ERROR=ERROR+ABS(UNEW(i,j)-UOLD(i,j))
END SUBROUTINE ERROR_CALCULATION

SUBROUTINE RR(UNEW,U0,R)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP),DIMENSION(N,N),INTENT(IN)::U0,UNEW
    REAL(DP),DIMENSION(N,N),INTENT(OUT)::R
R(i,j)=U0(i,j)-UNEW(i,j)
END SUBROUTINE RR

END PROGRAM ELLIPTIC_SOLVER