Add sparse matrix interface

[majosm]

WIP

cc @lukeolson

[majosm]

@inducer I added a loopy kernel fallback for eager (code), and it's working when the array being multiplied is 1D, but I'm running into some trouble with array strides for higher dimensions. I set up the input array like this:

    u1 = actx.np.where(actx.np.less_equal(x, (a+b)/2), 0*x + 1, 0*x)
    u2 = actx.np.where(actx.np.greater_equal(x, (a+b)/2), 0*x + 1, 0*x)
    # u = actx.freeze_thaw(u1)  # 1D
    u = actx.freeze_thaw(actx.np.stack([u1, u2]).T)  # 2D

The strides are no longer C-like after transposing. Lazy doesn't seem to care about this, but eager doesn't like it; it makes my 1D diffusion solver blow up instantly. If I add order="F" to the GlobalArg(...) calls for "array" and "out", then it works fine. Is there a way I can fix this in general?

Edit: I guess the strides of u are only non-C-like in the eager case.

[majosm]

FWIW, using a dense matrix + actx.einsum instead shows the same behavior.

[majosm]

Upon closer inspection, it looks like the kernels may actually be OK. The problem occurs later in the test driver, when I do this:

        u = actx.freeze_thaw(u + dt*compiled_rhs(t, u))

(here u has order='F' and dt*compiled_rhs(t, u) has order='C' in the eager case). They don't add together correctly in eager, which seems to be due to inducer/pyopencl#681. So I guess the solution is just "don't mix orders".