Over the last couple of months I've been occasionally working on a package for iterative solution of linear systems. We already have several direct solvers for dense systems based on matrix factorizations. With this addition we would have solvers for large sparse systems, for systems where the matrix doesn’t need to be explicitly formed, or for systems where the solution is not required to be accurate to machine precision (for example in inexact Newton methods for the optimize package).
There are certainly some corners that need to be polished but I hope that it has reached a presentable state. The code currently lives at:
https://github.com/vladimir-ch/iterativeIf it is accepted here, I will then submit a PR to the gonum repository to discuss the code itself.
First of all, the package name. I settled on “iterative” very early on and so far I didn’t come up with a better name. It reads quite nicely: the main interface is iterative.Method, the main function is iterative.LinearSolve.
The overall design is similar to the optimize package. At the low level there are various solver types that implement the Method interface. It uses the reverse communication mechanism in order to avoid depending on how the system matrix is represented. Operations like (transposed) matrix times vector are delegated to the caller. In this way we can provide the iterative algorithms without having any implementation for sparse matrices available. At a higher level there is a convenience function LinearSolve which implements the main iteration loop similarly to optimize.Local. The convenience comes at the “cost” of concretising how the (transposed) matrix times vector (and other) operations are done.
The implementation of the various methods loosely follows “The Templates”
http://www.netlib.org/templates/double/. Actually, I used the version in SciPy where the bugs present in the netlib version are fixed. For the start I would have been happy with providing just CG and GMRES but to better understand what operations are needed, I did also BiCG and BiCGSTAB. I plan to focus on Krylov methods, so methods like SOR don’t fit into the current design.
Regarding the API, I’d like to get some opinions mainly about how to represent the matrix operations (and the preconditioner). At the moment there is a struct MatrixOps with two fields MatVec func(dst, x []float64) and MatTransVec func(dst, x []float64). The other obvious option is to make it an interface which is probably better. I have other smaller doubts (e.g., details of the reverse communication protocol: who should compute what, what operations are necessary), and there are many TODOs sprinkled over the code but I hope these will be raised and addressed during the review.