xORBDB6/xUNBDB6: fix a constant

[christoph-conrads]

Description

The commits listed below updated xORBDB6/xUNBDB6 with the aim of improving numerical stability by avoiding the use of squared values. With the exception of SORBDB6 these commits did not properly up a magic constant called ALPHA in the code. SORBDB6 was initially updated correctly before commit a015b21 inserted an incorrect value again.

The incorrect magic constant may lead to early returns with only partially orthogonalized output vectors.

Relevant commits:

• 54b3964
• 6479e0f
• 94419e8

Checklist

• The documentation has been updated.
• If the PR solves a specific issue, it is set to be closed on merge.

[christoph-conrads]

Input triggering the bug fixed in this PR can be found in the linked post: #917 (comment).

[langou]

Thanks @christoph-conrads. ALPHA = 1 / 100 is certainly safer than ALPHA = 1/10.

[christoph-conrads]

Thanks @christoph-conrads. ALPHA = 1 / 100 is certainly safer than ALPHA = 1/10.

ALPHA is now one tenth, that is ALPHA is larger. The constant is used in two places:

• to check if an early return is possible (see line 261ff. in CUNBDB6)
• to check if the vector x is practically zero after the second Gram-Schmidt iteration (see line 307ff. in CUNBDB6)

Moreover the value was not really changed. In the past it was used in comparisons involving squares of norms (e.g., as in ALPHASQ = 0.01; ALPHASQ ≤ ∥x∥²); the current code avoids the use of squares (e.g., as in ALPHA = 0.1; ALPHA ≤ ∥x∥).

[langou]

lapack/SRC/dorbdb6.f

Lines 255 to 263 in 8ecaaf9

	*
	* If projection is sufficiently large in norm, then stop.
	* If projection is zero, then stop.
	* Otherwise, project again.
	*
	IF( NORM_NEW .GE. ALPHA * NORM ) THEN
	RETURN
	END IF
	*

Thanks for correcting me.

ALPHA = 1/100 leads to a loss of accuracy of about 100x. ALPHA = 1/10 is better. Sure. Good. We can continue stir ALPHA up. A standard value, in the Gram-Schmidt context, is ALPHA = 1/sqrt(2). (But we are not in the Gram-Schmidt context.)

Since, this projection (maybe with reorthogonalization) is one and done, as opposed to repeatedly project columns one after the other as in the case in Gram-Schmidt, the value of ALPHA is a little different.

Taking an ALPHA small such as 1/10 should not be too much of an issue in this case. (Only 10x loss of accuracy at most.)

All in all though, I am not against setting ALPHA to ONE so as to force reorthogonalization in all cases.