# returns -∇² (discrete Laplacian, real-symmetric positive-definite) on n₁×n₂ grid
function ∇²(n₁,n₂)
    o₁ = ones(n₁)
    ∂₁ = spdiagm((-o₁,o₁),(-1,0),n₁+1,n₁)
    o₂ = ones(n₂)
    ∂₂ = spdiagm((-o₂,o₂),(-1,0),n₂+1,n₂)
    return kron(speye(n₂), ∂₁'*∂₁) + kron(∂₂'*∂₂, speye(n₁))
end

∇²(4,3)

full(∇²(4,3))

using PyPlot

spy(∇²(10,5))

A = ∇²(600,200) # a 120,000×120,000 matrix
b = rand(600*200)
@time A \ b
@time cholfact(A) \ b
@time lufact(A) \ b;

spy(sparse(cholfact(∇²(10,5))))
title("sparse Cholesky factor, AMD ordering")

if VERSION < v"0.4.0-dev+4289"
    error("You need a newer version of Julia to access user-specified sparse-Cholesky permutations")
end

spy(sparse(cholfact(∇²(10,5), perm=1:50)))
title("sparse Cholesky factor, natural ordering")

N = 10:10:400
m = N.^2
loglog(m, [nnz(cholfact(∇²(n,n), perm=1:n*n)) for n in N], "r*")
loglog(m, [nnz(cholfact(∇²(n,n))) for n in N], "bo")
loglog(m, m.^(3/2), "r-")
loglog(m, m .* log2(m), "b-")
xlabel(L"matrix size $m$")
ylabel("# nonzeros in Cholesky")
legend(["natural", "AMD", L"m^{3/2}", L"m \log_2(m)"], loc="upper left")

# return nested-dissection permutation for an m × n grid (in column-major order).  ndissect is an optional
# number of dissection steps before reverting to the normal ordering, and lm is an optional leading dimension
# (used for indexing computations).
function nest(m,n, ndissect=30, lm=m)
    if ndissect <= 0 || m*n < 5
        return reshape([i + (j-1)*lm for i=1:m, j=1:n], m*n)
    elseif m >= n
        msep = div(m,2)
        N1 = nest(msep-1,n, ndissect-1, lm)
        N2 = nest(m-msep,n, ndissect-1, lm) + msep
        Ns = msep + (0:n-1)*lm
        return [N1; N2; Ns]
    else
        nsep = div(n,2)
        N1 = nest(m,nsep-1, ndissect-1, lm)
        N2 = nest(m,n-nsep, ndissect-1, lm) + nsep*lm
        Ns = (1:m) + (nsep-1)*lm
        return [N1; N2; Ns]
    end
end

sort(nest(30,60)) == [1:30*60;]

spy(sparse(cholfact(∇²(10,5), perm=nest(10,5))))
title("sparse Cholesky factor, nested-dissection ordering")

m = 100
n = 120
nnz(cholfact(∇²(m,n), perm=1:m*n)), nnz(cholfact(∇²(m,n), perm=nest(m,n))),  nnz(cholfact(∇²(m,n)))

figure(figsize=(10,10))
n1, n2 = 10, 5
m = n1*n2

subplot(2,2,1)
F = cholfact(∇²(n1,n2), perm=1:m)
spy(sparse(F))
title("natural, $(nnz(F)/m^2*100)% fill", y=1.1)

subplot(2,2,2)
F = cholfact(∇²(n1,n2), perm=nest(n1,n2, 1))
spy(sparse(F))
title("1 dissection, $(nnz(F)/m^2*100)% fill", y=1.1)

subplot(2,2,3)
F = cholfact(∇²(n1,n2), perm=nest(n1,n2, 2))
spy(sparse(F))
title("2 dissections, $(nnz(F)/m^2*100)% fill", y=1.1)

subplot(2,2,4)
F = cholfact(∇²(n1,n2), perm=nest(n1,n2, 3))
spy(sparse(F))
title("3 dissections, $(nnz(F)/m^2*100)% fill", y=1.1)