Issue 4136: Specify behavior of [linalg] Hermitian algorithms on diagonal with nonzero imaginary part

4136. Specify behavior of [linalg] Hermitian algorithms on diagonal with nonzero imaginary part

Section: 29.9.3 [linalg.general] Status: New Submitter: Mark Hoemmen Opened: 2024-08-09 Last modified: 2024-08-10

Priority: Not Prioritized

View all issues with New status.

Discussion:

Mathematically, the diagonal elements of a Hermitian matrix must all have zero imaginary part (if they are complex numbers). 29.9.3 [linalg.general] paragraphs 3 and 4 govern the behavior of [linalg] functions that operate on "symmetric," "Hermitian," or "triangular" matrices. All such functions only access the specified triangle (lower or upper) of the matrix. Whatever is in the other triangle of the matrix doesn't matter; it's not even accessed. That gives well-defined behavior for "symmetric" and "triangular" matrices. However, both triangles include the diagonal. What should the "Hermitian" functions do if they encounter a diagonal element with nonzero imaginary part?

The current wording says that both the operation performed and the matrix itself are Hermitian, but does not clarify what happens if the latter is not true. For example, 29.9.14.3 [linalg.algs.blas2.hemv] says that hermitian_matrix_vector_product performs a

Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A (29.9.3 [linalg.general]).

Language like this appears in the specifications of all the functions whose names start with hermitian. The implication is that if the diagonal has an element with nonzero imaginary part, then the matrix is not Hermitian and therefore a precondition of the function has been violated. The result is undefined behavior.

We can get rid of this undefined behavior by defining what happens in this case. It turns out that Chapter 2 of the BLAS Standard already does this: It says that the Hermitian algorithms do not access the imaginary parts of diagonal elements. The reference Fortran BLAS implementations of CHEMV (single-precision Complex HErmitian Matrix-Vector product) and ZHEMV (double-precision complex HErmitian Matrix-Vector product) follow the BLAS Standard. CHEMV uses real(A(j,j)) and ZHEMV uses dble(A(j,j)), which means that the BLAS routines only access the real part of each diagonal element.

The clear design intent of P1673R13 was to imitate the behavior of the BLAS Standard and reference BLAS unless otherwise specified. Thus, we propose to specify this behavior in the wording; we do not think this requires LEWG re-review.

In my view, it's fine to retain the existing wording like that in 29.9.14.3 [linalg.algs.blas2.hemv], that refers to a "Hermitian matrix". The matrix defined by the revised [linalg.general] instructions is unambiguously a Hermitian matrix. Just like the "other triangle" of a symmetric or triangular matrix does not affect the behavior of [linalg] symmetric resp. triangular algorithms, the wording fix here will ensure that any imaginary parts of diagonal elements will not affect the behavior of [linalg] Hermitian algorithms.

Proposed resolution:

This wording is relative to N4988.

Modify 29.9.3 [linalg.general] as indicated:

-4- For any function F that takes a parameter named t, t applies to accesses done through the parameter preceding t in the parameter list of F. Let m be such an access-modified function parameter. F will only access the triangle of m specified by t. For accesses of diagonal elements m[i, i], F will use the value real-if-needed(m[i, i]) if the name of F starts with hermitian. For accesses m[i, j] outside the triangle specified by t, F will use the value […]