# L-29 MCS 572 Fri 1 Nov 2024 : jacobi.jl

# Defines the method of Jacobi.

using LinearAlgebra
using Printf

"""
    function jacobi(A::Matrix{Float64},
                    b::Vector{Float64}, x0::Vector{Float64},
                    N::Int, tol::Float64, verbose=true)

does at most N iterations of Jacobi matrix, to solve A*x = b,
starting at the vector x0.  Stops when the norm of the update
vector is less than the tolerance tol.
If verbose, then the norm of each update is written to screen.
"""
function jacobi(A::Matrix{Float64},
                b::Vector{Float64}, x0::Vector{Float64},
                N::Int, tol::Float64, verbose=true)
    nit = N
    x = deepcopy(x0)
    D = Diagonal(A)
    for k=1:N
        dx = D^(-1)*(b - A*x)
        x = x + dx
        ndx = norm(dx)
        if verbose
            sdx = @sprintf("%.2e", ndx)
            println(sdx)
        end
        if ndx < tol
            return (k, x)
        end
    end
    return (nit, x)
end

"""
    function main()

tests Jacobi's method on a diagonally dominant matrix.
The dimension is taken from the command line arguments if provided.
"""
function main()
    if length(ARGS) > 0
        dim = parse(Int, ARGS[1])
    else
        dim = 10
    end
    I = ones(dim, dim)
    A = I + dim*Diagonal(I)
    b = 2*dim*ones(dim)
    println(A)
    x = A\b
    println(x)
    x0 = zeros(dim)
    N = dim*dim
    tol = 1.0e-4
    n, z = jacobi(A, b, x0, N, tol)
    println(n)
    println(z)
end

main()