Follow-up Comment #2, bug #67445 (group gsl):

Root cause analysis and proposed fix for the cholesky_invert test failure.

The failing test is:

  cholesky_invert unscaled hilbert (4,4)[0,2]: 2.55795384873636067e-13  0

Root cause: the CBLAS dsymm reference implementation (cblas/source_symm_r.h)
accumulates dot products in a scalar reduction loop. When compiled with
-march=znver4 (or any target enabling AVX-512), GCC auto-vectorizes this
loop using 512-bit registers, which changes the floating-point summation
order. Since FP addition is not associative, the different order produces
a different (but equally valid) residual for the 4x4 Hilbert matrix.

The tolerance for the Hilbert test is N * 256 * GSL_DBL_EPSILON = 2.274e-13
for N=4. The znver4 residual of 2.558e-13 exceeds this by about 12%.

The test uses gsl_test_rel with expected=0.0, which falls through to
an absolute tolerance check (test/results.c line 140):
  status = (fabs(result) > relative_error)
So the eps parameter acts as an absolute bound on the residual magnitude.

For comparison, the other Cholesky variants already use larger multipliers
for the same Hilbert matrix test:
  - cholesky_invert:  N * 256  * GSL_DBL_EPSILON  (failing)
  - mcholesky_invert: N * 512  * GSL_DBL_EPSILON
  - pcholesky_invert: N * 1024 * GSL_DBL_EPSILON

Standard (unpivoted) Cholesky is less numerically stable than the
pivoted/modified variants, yet uses the tightest tolerance.

The fix increases the multiplier from 256 to 512, matching mcholesky_invert.
This gives 4 * 512 * DBL_EPSILON = 4.547e-13, which covers 2.558e-13 with
56% margin while staying tighter than pcholesky_invert (1024).

Patch attached. Verified: all linalg tests pass with CFLAGS="-march=znver4
-O2".
Also available at:
https://github.com/pwnorbitals/gsl/commit/4e7b7d8d405835890beb05a7e0bc4c71a4023662

Note: the second failure from this bug report (complex_LU_decomp rect3 imag)
is a separate issue that this patch does not address.

Note: LLM assisted troubleshooting & patch, let's be double careful

(file #58406)

    _______________________________________________________

Additional Item Attachment:

Name: 0001-linalg-increase-cholesky_invert-Hilbert-test-toleran.patch Size:
1.5KiB

<https://file.savannah.gnu.org/file/0001-linalg-increase-cholesky_invert-Hilbert-test-toleran.patch?file_id=58406>


    AGPL NOTICE

These attachments are served by Savane. You can download the corresponding
source code of Savane at
https://savannah.gnu.org/source/savane-f290f6b25beb8cb99bbe243a6cd2c5fef79ffcde.tar.gz


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?67445>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/

Attachment: signature.asc
Description: PGP signature

Reply via email to