search

kokkos / stdBLAS

Reference Implementation for stdBLAS
Other
118 stars 22 forks source link

Fix symmetrical matrix vector product #278

Closed iburyl closed 1 month ago

iburyl commented 1 month ago

Algorithms for symmetric_matrix_vector_product and hermitian_matrix_vector_product double counted the diagonal. This fix is to count diagonal only once.

dalg24 commented 1 month ago

Mark do you want to merge as is or try to come up with some test that would expose that bug?

iburyl commented 1 month ago

Wait-wait-wait, do not push this PR just yet!

Found one more inconsistency with BLAS implementation, which is not specifically documented in P1673. Thus the correct implementation depends on the answer.

BLAS hemv ignores imaginary part of the A diagonal elements, considering them purely real. https://www.netlib.org/lapack/explore-html/db/d17/group__hemv_ga8e55e480b23945c09c35e09935fce058.html#ga8e55e480b23945c09c35e09935fce058 See lines with temp2 = zero

This is not the case in this PR.

What is the assumed behavior in case of P1673 for the imaginary part of diagonal elements?

iburyl commented 1 month ago

I assume that BLAS behavior is the correct one, because Hermitian matrices should have Real diagonal by definition. Silently ignoring part of the diagonal is something not necessarily in spirit of std C++, but I have no idea what we can do otherwise... Thus, one more commit above.

mhoemmen commented 1 month ago

@iburyl wrote:

What is the assumed behavior in case of P1673 for the imaginary part of diagonal elements?

That's effectively UB. [linalg.algs.blas2.hemv] 5 says: "These functions perform an overwriting Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A." If A's diagonal has elements with nonzero imaginary parts, then A is not Hermitian.

Chapter 2 of the BLAS Standard ( https://www.netlib.org/blas/blast-forum/chapter2.pdf ) says: "The imaginary part of the diagonal entries of the matrix operand are supposed to be zero and should not be referenced." The various error checking specifications of the BLAS Standard say nothing about what happens if they are nonzero. In the reference BLAS, ZHEMV invokes dble on the diagonal element, which is like taking the real part.

https://www.netlib.org/lapack/explore-html/db/d17/group__hemv_ga8e55e480b23945c09c35e09935fce058.html#ga8e55e480b23945c09c35e09935fce058

I would interpret "should not be referenced" as "it's UB to read those values." The reference BLAS appears to follow the same interpretation.

@dalg24 wrote:

Mark do you want to merge as is or try to come up with some test that would expose that bug?

Thanks for asking! It would be excellent if Ilya has the time to add a test, but if he doesn't, I'd rather just merge the fix first. It's clearly broken. We talked about this yesterday, and I don't have time to add a text over the next few days.

iburyl commented 1 month ago

I'd rather not count on myself in the short term regarding the test.