[SciPy-dev] implement reorthogonalization method in gmres

[Pauli Virtanen]

to, 2010-02-25 kello 14:16 -0500, Darcoux Christine kirjoitti:
> Thank you for your quick answer with a pointer to the litterature.
> 
> The result in the paper you cite assumes that the matrix is
> well-conditioned, and I'm not sure I can count on that in a
> general-purpose code which many people will use.

Yes, it's a good point that the situation might not be the same for near
singular matrices. Perhaps someone looked into that...

> The books "Solving nonlinear equations with Newton's method" and
> "Iterative methods for linear and nonlinear equations" (some parts are
> availlable from google book) by C. T. Kelley gives a good presentation
> on the importance of orthogonality. It also show a simple 3x3 example
> where paying attention to orthogonality is very important.

I note that the book is from ca 1995, whereas the papers are newer. For
the particular example in the book, 

        [ 0.001  0  0;  0  0.0011  0; 0 0 1e4 ]

MGS (which we currently have) required an additional step at the final
stage of the convergence, as compared to GS + selective
reorthogonalization. Up to the last step, the convergence histories were
identical.

> Checking for orthogonality and reorthogonalizing the way there author
> of these books do it costs almost nothing (see the link to the code in
> my previous post). Matlab's "official" GMRES code uses Householder
> reflections to build the basis and gets orthogonality automatically,
> but at double the cost of Gram-Schmidt. The Sandia Trilinos code uses
> classical Gram-Schmidt for good parallel performance and
> reorthogonalizes at every iteration, so also costs double.
> The way proposed by Kelley costs nothing because the
> reorthogonalization happens rarely, so I do not see any reason not to
> do it. Of course, this could be added as an option so users will
> decide if they want to use it.

Well, as I said, if someone implements it + writes tests, it will get
in. If you are not going to implement it, please submit an enhancement
request at

        http://projects.scipy.org/scipy

to keep track of the idea. I don't have the time to work at this right
now, though, but later is possible.

Out of interest: do you have an actual matrix at hand where you know
that MGS + selective reorthogonalization will give significantly better
results than plain MGS, or a reference that shows an example like that?

-- 
Pauli Virtanen