kokkos / stdBLAS

Reference Implementation for stdBLAS
Other
127 stars 22 forks source link

"symmetric_matrix_rank_1_update" can only be used for positive(semi) definite updates #262

Closed yigiter closed 9 months ago

yigiter commented 10 months ago

(This issue is probably a duplicate of #138 / #136. Yet, I still wanted to raise the following issue in case completely removing "alpha" has been decided as the final solution for #138)

BLAS definition for Hermitian/symmetric rank-1 update (xSYR) is as follows (alpha is real number which may be less than 0): A := alpha*x*x**T + A

However, symmetric_matrix_rank_1_update(x,A,t) in stdBlas implements this as follows: A := x*x**T + A

I believe the intention is that the users should define the scale (alpha) using the scaled() function. For instance (assuming alpha=2): symmetric_matrix_rank_1_update(scaled(sqrt(2.0),x),A,t)

However, with such a usage, it is not possible to utilize a negative definite symmetric rank-1 update. (For instance, how can symmetric_matrix_rank_1_update be used when alpha=-2?)

For general rank 1 updates, this may not be an issue as the scaled() function can be used as intended on a single argument. However, for the symmetric and hermitian rank-1 updates, an explicit alpha argument seems to be necessary in order to be compatible with xSYR.

mhoemmen commented 9 months ago

@yigiter Thanks for the bug report!

P1673R13 (the final approved version) defines overloads with and without alpha for the reason you mentioned. The reference implementation needs a few updates to bring it in line with. We've had these overloads in R12 but it's possible that the reference implementation isn't correctly constrained or is missing the overloads. We'll take a look as soon as we get the chance. Thanks!

mhoemmen commented 9 months ago

This is related to #261 , in that the ExecutionPolicy overloads might not be constrained enough for the compiler to resolve the overloads that don't take vs. take alpha.

mhoemmen commented 9 months ago

@yigiter My PR #263 (merged today) ended up fixing this issue as well. I added a test for symmetric_matrix_rank_1_update that calls it both with and without an alpha. Thanks for contributing this bug report!