module Blocksys

export load_matrix, load_vector, compute_vector, save_vector, gaussian_elimination!

bcomputed = true

function load_matrix(filename::String)
    open(filename) do f
        line = readline(f)
        (n, l) = split(line)
        (n, l) = (parse(Int64, n), parse(Int64, l))
        Is = Array{Int64, 1}()
        Js = Array{Int64, 1}()
        Vs = Array{Float64, 1}()
        println("n: ", n, ", l: ", l)
        while !eof(f)
            line = split(readline(f))
            push!(Is, parse(Int64, line[1]))
            push!(Js, parse(Int64, line[2]))
            push!(Vs, parse(Float64, line[3]))
        end
        return (sparse(Is, Js, Vs), n, l)
    end
end

function load_vector(filename::String)
    open(filename) do f
        line = readline(f)
        n = parse(Int64, line)
        b = Vector{Float64}()
        while !eof(f)
            line = readline(f)
            push!(b, parse(Float64, line))
        end
        global bcomputed = false
        return b
    end
end

function gaussian_elimination!(A::SparseMatrixCSC{Float64, Int64}, n::Int64, l::Int64)
    for k = 1:n-l
        println("col: ", k)
        rowToReset = k + l +l - 1
        for i = k+1:rowToReset
            A[i, k] = 0
        end
    end
    show(IOContext(STDOUT), "text/plain", full(A))
end

function compute_vector(A::SparseMatrixCSC{Float64, Int64}, n::Int64)
    b = zeros(Float64, n)
    vals = nonzeros(A)
    rows = rowvals(A)
    for i = 1:n
        for j in nzrange(A, i)
            b[rows[j]] += vals[j]
        end
    end
    global bcomputed = true
    return b
end

function save_vector(x::Vector{Float64}, filename::String)
    open(filename, "w") do f
        length = size(x, 1)
        if bcomputed
            exactX = ones(length)
            relativeError = norm(exactX - x) / norm(exactX)
            println(f, relativeError)
        end
        for i = 1:length
            println(f, x[i])
        end
    end
end


loaded = load_matrix("A.txt")
println("n: ", loaded[2], ", l: ", loaded[3])
show(IOContext(STDOUT), "text/plain", full(loaded[1]))
println()
gaussian_elimination!(loaded[1], loaded[2], loaded[3])
#=
tic()
bvector = compute_vector(loaded[1], loaded[2])
toc()
println("b: ", bvector)
=#

loadedB = load_vector("b50k.txt")
#save_vector(loadedB, "test.txt")
end