Simplify construction of 3D Poisson matrix

Author: haampie

test/cg.jl

@@ -4,7 +4,7 @@ using LinearMaps
 
 srand(1234321)
 
-include("getDivGrad.jl")
+include("poisson_matrix.jl")
 
 facts("cg") do
 
@@ -39,7 +39,7 @@ context("Small full system") do
 end
 
 context("Sparse Laplacian") do
-    A = getDivGrad(32,32,32)
+    A = poisson_matrix(Float64, 32, 3)
     L = tril(A)
     D = diag(A)
     U = triu(A)

test/getDivGrad.jl

@@ -0,0 +1,12 @@
+function poisson_matrix{T}(::Type{T}, n, dims)
+    D = second_order_central_diff(T, n);
+    A = copy(D);
+
+    for idx = 2 : dims
+        A = kron(A, speye(n)) + kron(speye(size(A, 1)), D);
+    end
+
+    A
+end
+
+second_order_central_diff{T}(::Type{T}, dim) = convert(SparseMatrixCSC{T, Int}, SymTridiagonal(fill(2 * one(T), dim), fill(-one(T), dim - 1)))