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.