A dumb thought: technically scipy has a “solve_banded” function that does a banded solve. He could easily recast his problem into a single big diagonal problem, I guess, just with some silly explicit zeros added. I wonder how the performance of that would compare, to iterating over a bunch of tiny solves.
Of course, it would be nice if scipy had a block diagonal solver (maybe it does?). But yea, I mean, of course it would be nice if my problem was always built in functionality of the library, lol.
Maybe a bsr_matrix could work.