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