! Filename : SOR.f
! Edited by SHIN, 2004. 5. 3
! To solve A*x = b using SOR Method
 
! -----------------------------------------------------------
Program Main
  Implicit None
  real, Allocatable :: A(:,:), b(:), x(:), ex(:)
  Integer :: i,j,n,itr=10,maxit=3000
  Real :: omega,l2er,h1er,tol=.1**10
  Real :: Start_t, end_t
  Real :: Dot_Product, Matmul, l2norm, Matnorm
! -----------------------------------------------------------
  
  Print *, ' Input the dimension of Problem n = '
  Read *, n
  if (mod(n,2)==0) n=n+1

  Allocate( A(n,n),b(n),x(n),ex(n))
! -- Give an example of A, b, x0 --
  A = 0.0 ; x = 0.0 ; 
  Do i=1,n
    A(i,i-1)=-1.0; A(i,i)=2.0; A(i,i+1)=-1.0
  Enddo
!  Call Random_Number(ex); ex = 90*ex+10 ! Random Exact solution
  ex = 2.0
  b = Matmul(A,ex)                      ! Right Hand Side with ex
! -- Give an example of A, b, x0 --

  Print *, 'Input the relaxation parameter omega = '
  Read *, omega
  Call CPU_Time(start_t)
! -- Call SOR method
    print *, ' Iteration  l2-error   h1-error '
SorDo : Do j=1,maxit
    call SOR(A,b,x,n,omega,itr)
    l2er = l2norm(x-ex,n) ; h1er = Matnorm(A,x-ex,n)
    write(*,'(i5,2x,D14.6,2x,D14.6)') itr*j ,l2er ,h1er
    if (l2er < tol) Exit
  Enddo SorDo
  Call CPU_Time(end_t)
  print *,'=========================================================='
  Print *, "  CPU-Time  Dimension  Iterations     Error       omega"
  Write(*,"(F10.3, 4x, I5, 8x, I5, 2x, 1pQ14.5, 2x, 0pF6.2)") &
          & end_t-start_t, n, j*itr, l2er, omega
  print *,'=========================================================='

End Program Main

!-----------------------------------------------------------
Subroutine SOR(A,b,x,n,omega,itr)
  Implicit None 
  Integer :: n,i,itr,k
  Real :: omega, A(n,n), x(n), b(n)
!-----------------------------------------------------------
 Do k=1,itr
  Do i=1,n
    x(i)=(1-omega)*x(i)+(omega/A(i,i))*( &
               b(i)-A(i,i-1)*x(i-1)-A(i,i+1)*x(i+1) )
  Enddo 
 Enddo 
End Subroutine SOR  

!-----------------------------------------------------------
Real Function l2norm(x,n)
  Implicit None ; Integer n ; Real x(n)
!-----------------------------------------------------------
  l2norm = Dot_Product(x,x); l2norm = sqrt(l2norm)
End Function l2norm
!-----------------------------------------------------------
Real Function Matnorm(A,x,n)
  Implicit None ; Integer n ; Real A(n,n), x(n)
!-----------------------------------------------------------
  Matnorm = Dot_Product(x,Matmul(A,x)); Matnorm = sqrt(Matnorm)
End Function Matnorm