[petsc-dev] suggestion of change in line 93 for petsc 3.6.0 for ./src/ksp/ksp/impls/gmres/borthog2.c

[Barry Smith]

   Julien,

    Thanks for the complete report. I have fixed the constant in the maint, master, and next branches. You are also right about the constant I used, it came from my understanding of how ARPACK handled it.

   Barry

> On Jun 12, 2015, at 2:53 PM, Langou, Julien <Julien.Langou at ucdenver.edu> wrote:
> 
> Hi Barry, 
> 
> This is Julien from CU Denver. 
> 
> Not sure to whom to report this suggestion of change so I am pushing to you. 
> 
> Line 93 of ./src/ksp/ksp/impls/gmres/borthog2.c reads:
> 
> if (wnrm < 1.0286 * hnrm) {
> 
> I suggest to change it to:
> 
> if (wnrm < hnrm) {
> 
> The constant 1.0286 is kind of weird in this context. 
> 
> It is hard for me to explain. 
> 
> I have a few slides about this that I have attached. This is extracted from a talk on the numerics of the Gram-Schmidt algorithm. 
> 
> I think I can trace the 1.0286 in PETSc ( instead of 1 ) back to ARPACK. 
> 
> ARPACK intended to use SQRT( 2 ) — which is a standard value, that lots of numerical codes tend to agree upon. Someone in the ARPACK team hard-coded SQRT( 2 ) as 0.7170 which is unfortunate, since the correct value is closer to 0.7071. I think someone in PETSC copied ARPACK orthogonalization scheme, the computation is not done with the same quantities so what corresponds to 0.7170 is actually correctly computed as 1.0286, but the if one were to take the intended sqrt( 2 ) value, one should simply get 1 as opposed to 1.0286.
> 
> So I suggest to change the 1.0286 to a 1 in the PETSc code. 
> 
> For reference, we (Giraud and Langou, SISC, 2003) proved that no constant will really work. So if you take wnrm < 1/2 * hnrm, you can expect better orthogonality but we prove that 1 is not enough, 1/2 is not enough, 1/3 is not enough, etc. We also give a criteria that guarantees orthogonality. (up to numerical accuracy) but, to some extent, our criteria is most of the time too conservative and 1 (which is what you intended with 1.0286) is OK and I believe it is what people use in practice. This means that if the angle between the vector to be made orthogonal and the orthogonal set is less than 45 degrees then two orthogonalization steps will be made whereas, if the angle is more than 45 degrees, then only one orthogonalization will be made. That’s all good with me. 
> 
> Best wishes, 
> Julien.
> <Overview GS.pdf>