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 <[email protected]>
> 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>
