|
New iterative sparse solver: LGMRES
Implement a new pure-Python sparse solver LGMRES based on the article
A.H. Baker and E.R. Jessup and T. Manteuffel,
SIAM J. Matrix Anal. Appl. 26, 962 (2005).
http://dx.doi.org/10.1137/S0895479803422014
Preprint:
http://www.osti.gov/energycitations/product.biblio.jsp?osti_id=15020368
http://www.osti.gov/energycitations/servlets/purl/15020368-IonFj3/15020368.pdf
It is designed to avoid the "stagnation" problem in restarted GMRES, and so
often outperforms GMRES by some measure, or at least is not much worse. For
example (try it out yourself, scipy/sparse/linalg/isolve/tests/demo_lgmres.py):
MatrixMarket problem SPARSKIT/drivcav/e05r0200
Invert 236 x 236 matrix; nnz = 5856
GMRES(100): 12372 matvecs, residual 1.88541625334e-14
LGMRES(88,6) [same memory req.]: 8189 matvecs, residual: 8.78114597227e-15
LGMRES(94,6) [same subspace size]: 7031 matvecs, residual: 6.26488040431e-15
Another advantage in this algorithm is that you can supply it with "guess"
vectors that augment the Krylov subspace; if the solution lies close to the
span of these vectors, the algorithm converges faster. This can be useful if
several very similar matrices need to be inverted one after another, such as in
Newton-Krylov iteration where the Jacobian matrix often changes little in the
nonlinear steps.
Profiling shows that even though the least-squares problem is solved with
`lstsq` and not on-line with Givens rotations, this accounts to only ~ 35% of
the total time even in the above problem with small NNZ (and large subspace
size M=100). So I believe optimizing this part is of low priority and can be
left for later
|