# L-11 MCS 471 Fri 16 Sep 2022 : sor.jl
# Illustrates the method of successive over-relaxation

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

using LinearAlgebra

"""
    sor(mat::Array{Float64},rhs::Array{Float64},
        sol::Array{Float64},wgt::Float64=1.25,
        maxit::Int=100,tol::Float64=1.0e-8) 

Runs successive over-relaxation on the linear system with coefficient matrix
in mat, with right hand side vector in rhs, and a start solution in sol,
with weight wgt.
Running stops if the maximum number of iterations in maxit is reached,
or if the norm of the correction is less than the given tolerance.

Returns (solution, numit, nrmdx, fail), the computed solution,
the number of iterations numit, an estimate for the forward error nrmdx,
and fail is true if the given tolerance was not reached.
"""
function sor(mat::Array{Float64},rhs::Array{Float64},
             sol::Array{Float64},wgt::Float64=1.1,
             maxit::Int=100,tol::Float64=1.0e-8) 
    nbrows, nbcols = size(mat)
    result = deepcopy(sol)
    deltax = zeros(nbrows, 1)
    numit = 0; nrmdx = 1;
    while numit < maxit
        numit = numit + 1
        for i=1:nbrows
            deltax[i] = wgt*rhs[i]
            for j=1:nbcols
                deltax[i] = deltax[i] - mat[i,j]*wgt*result[j]
            end
            deltax[i] = deltax[i]/mat[i,i]
            result[i] = result[i] + deltax[i]
        end
        nrmdx = norm(deltax)
        strdx = @sprintf("%.2e", nrmdx)
        println("||dx|| = $strdx")
        if norm(deltax) <= tol
            return (result, numit, nrmdx, false)
        end
    end
    return (result, numit, nrmdx, true)
end