Fix `find()` const-cast bug and diagmat bug

I'm the master of saying one thing and then immediately doing another. I thought I was going to go to bed, and here we are instead. This handles #30 (closed) and another bug I found earlier today.

• find() and variants would fail to call compute_n_rows() first in order to compute the size of a matrix. A simple const_cast fixes this, and I added tests for it.

• It turns out the other place I encountered #30 (closed) was the exact same thing, but for op_diagmat instead of op_diagmat2. The solution is simply to remove the overload in fn_diagmat.hpp, then add special handling for the operation in op_diagmat2::apply().

I tested on the CUDA backend only; I'd like to see it pass tests (or at least not fail tests relevant to this MR) before merge. Then a tag for release.