Very interesting, then why does the following code work?
using PyCall
@pyimport pyamg
@pyimport scipy.sparse as scipy_sparse
# generate 2D laplacian
N = 10
L1 = spdiagm((-ones(N-1), 2*ones(N), -ones(N-1)), (-1,0,1), N, N) * N^2
B = kron(speye(N), L1) + kron(L1, speye(N))
# load into Python
B_py_csr = scipy_sparse.csr_matrix(B)
# create multi-grid solver
ml = pyamg.ruge_stuben_solver(B_py_csr)
# solve with random RHS
b = ones(size(B,1))
x_py = ml[:solve](b, tol=1e-10)
# check result
println("|x - x_py| = ", norm(x_py - (B \ b), Inf))
After your comment I looked more closely. The line `B_py_csr =
scipy_sparse.csr_matrix(B)` takes **a lot** of time when N is large. Is it
possible that it converts B into a full matrix, then loads it into Python,
then generates the sparse matrix from that?
Christoph
