reordering

[gdmcbain]

In reviewing the example of forced convection #162, some suggestions were made to reduce the five-minute run-time. These concern the Navier–Stokes solver rather than forced convection per se, so new issues are being launched to take them up. The first suggestion was:

Maybe precalculating Reverse Cuthill-McKee reordering for the rows and columns of self.S would help a bit but can't say for sure.

Does reordering help in general? Does it help specifically with meshes imported from Gmsh and systems solved with scipy.sparse.linalg.solve? Does it depend on the block-structure of the matrix? (The Stokes problem as assembled in ex24 has the classic symmetric semidefinite saddle-point structure of the Taylor–Hood velocity–pressure mixed finite elements.) What about iterative solvers like GMRES as in the second suggestion?

The reordering is implemented in skfem.utils.rcm, but is not yet used in the examples.

[gdmcbain]

I hadn't bothered much about reordering hitherto as I recall a while ago reading a reply from one of the Gmsh authors on why Gmsh didn't reorder nodes. The rationale was that the way to do this was so dependent on the solver that it was best left to the latter. I got the impression that many modern canned solvers knew what was best for them in terms of order and took care of it themselves. I haven't looked into this for scipy.sparse.linalg.solve or the related Krylov solvers.

[gdmcbain]

The quote was (C. Geuzaine, 2017-07-28 [Gmsh] Mesh ordering)

No, we leave this to each solver, as each has different requirements.

[kinnala]

I think you're correct. If you plan to reuse a single LU-decomposition then you should reorder because it leads to less fill-in -> less memory usage and faster back substitution. But here you have a different matrix each time so it does not help.

[gdmcbain]

Ah. Although the matrix changes each time (each Reynolds number and each Newton iteration within that), I had been thinking that the sparsity would be the same so that the reordering could be shared, or even applied to the mesh (and then the basis regenerated). I'll leave this for the moment anyway.

[kinnala]

What you say is most certainly true, but I suppose scipy.sparse.linalg.solve will compute some (possibly better) reordering quite quickly before each solve anyways.

[gdmcbain]

The issue of reordering arose in the discussion of ‘Laplace “without” boundary conditions’ #592; in comparing ostensibly the same matrix with that generated in another code, the other had a tidier sparsity pattern.

[kinnala]

I suppose it could be possible to have an option to reorder mesh during initialization as an attempt to reduce fill in. However, I think it should not be enabled by default because it will cause surprising interpretation of the rows/cols of the resulting matrices.