Description

xORBDB5 and xORBDB6 compute vector norms differently causing disagreement on when a vector is numerically zero.

Given an isometric matrix Q and a vector x, xORBDB5 computes a vector x' that is orthogonal to the span of the columns of Q. (Q isometric means that the following two properties hold: Q has at most as many columns as rows and Q^* Q = I). Internally xORBDB5 calls xORBDB6. Given an isometric matrix Q and a vector x, xORBDB6 projects x onto the complement of the column span of Q. The code contains only the ominous comment that

Kahan's "twice is enough" criterion

is applied. Probably this means that xORBDB6 uses at most two iterations of classical Gram-Schmidt orthogonalization to compute its results. This approach is known as CGS2, see BarlowS2011 or GiraudLR2002 and it matches the matrix-vector multiplications inside this function. The problem are the computation of the vector norm.

For a xORBDB5, the norm of the vector computed by xORBDB6 is computed with SNRM2:

For xORBDB6, the vector norm is computed as follows:

Consider the input to xORBDB6 below:

c := 2^-45
μ := c ε

    | 0  1 -μ |
Q = | 0  0 +μ |
    | 1  0  0 |
    | 0  μ  1 |

With x = e_2 (i.e., the second standard basis vector), the xORBDB6 variable NORMSQ2 is zero whereas SNRM2 is nonzero causing xORBDB5 to return an incorrect vector.

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue