Try out block Jacobi iteration

[sampotter]

I'm pretty sure a good solver for the radiosity system will be a "block Jacobi iteration". Going to try to sketch out what this in the issue before trying to implement it.

For a system Ax = b, the Jacobi iteration separates A = D + L + U, D = diag(A), splitting A into D and L + U terms, then relaxing. This works very reliably for this system (seems to converge to machine precision in at most 8 iterations). It's also pretty simple since D = rho*I for the radiosity system.

A "block Jacobi iteration" would separate A = D + L + U where now D = diag(D1, ..., Dn), where the Dis are diagonal blocks at a certain depth. The idea here would be to relax as before, i.e. x <- D\(b - (L + U)x), except now D\ can be done by recursively applying the block Jacobi iteration, since each Di is a smaller radiosity system. At a fine enough level we can just use a sparse direct solver and avoid iterating. This leads to something reminiscent of a multigrid V-cycle (maybe there's another name for this in this context...). It seems likely that we would only need one or two of these "V-cycles" to converge fully.