PROGRAM SOR
    implicit none
    ! example
    double precision :: A(4,4) = reshape( (/ &
        1., -1., 0., 0., &
        -1., 2., -1., 0., &
        0., -1., 2., -1., &
        0., 0., -1., 2. &
        /), shape(A) )
    double precision :: b(4) = (/ 1., 0., 0., 0. /)
    double precision :: x(4) = 0
    double precision :: x_before(4) = 0
    double precision :: delta(4) = 0
    double precision :: threshold = 0.000000000001

    double precision w, temp_r
    integer i, j, dim
    
    w = 1.6
    dim = 4

    do while( .true. )
        do i = 1, dim
            temp_r = 0
            temp_r = temp_r + b(i)
            do j = 1, i - 1
              temp_r = temp_r - A(j, i) * x(j)
            end do
            do j = i + 1, dim
              temp_r = temp_r - A(j, i) * x_before(j)
            end do
            temp_r = temp_r/A(i,i)

            x(i) = x(i) + w * (temp_r - x(i))
            delta(i) = x_before(i) - x(i)
            x_before(i) = x(i)
        end do

        temp_r = 0
        do i = 1, dim
            temp_r = temp_r + delta(i) * delta(i)
        end do
        temp_r = sqrt(temp_r)
        if (temp_r < threshold) then
            exit
        end if   
    end do

    print *, '################'
    print *, 'answer'
    do i = 1, dim
        print *, x(i)
    end do

END PROGRAM