function getDivGrad(n1,n2,n3)
    # the Divergence
    D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
    D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
    D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))
    # DIV from faces to cell-centers
    Div = [D1 D2 D3]

    return Div*Div';
    end

function ddx(n)
# generate 1D derivatives
   return d = spdiagm((-ones(n),ones(n)),[0,1],n,n+1)
end

K = getDivGrad(60,60,60); b = rand(size(K,1)); c = similar(b);

@time for i in 1:1000 A_mul_B!(c,K,b) end;