Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add spr! to LinearAlgebra.BLAS #42830

Merged
merged 3 commits into from
Nov 22, 2021
Merged

add spr! to LinearAlgebra.BLAS #42830

merged 3 commits into from
Nov 22, 2021

Conversation

bjarthur
Copy link
Contributor

i used spmv! as a template

@ViralBShah ViralBShah added the linear algebra Linear algebra label Oct 27, 2021
NEWS.md Outdated Show resolved Hide resolved
@bjarthur
Copy link
Contributor Author

bjarthur commented Oct 28, 2021

PackedArrays.jl btw is an as of now private package which i plan to register and make public soon hopefully which implements a new SymmetricPacked type for julia. i prototyped the spr! wrapper there and inadvertently forgot to change the tests when copying into LinearAlgebra.

@bjarthur
Copy link
Contributor Author

i don't think any of the failing tests are this PRs fault

Comment on lines +1181 to +1182
α::Real, x::Union{DenseArray{T}, AbstractVector{T}},
AP::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasReal}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a usecase for x::DenseArray or would it make more sense to restrict it to just AbstractVector?

For AP, on the other hand, I don’t think we should allow AbstractVector without checking that its stride is compatible, i.e., chkstride1(AP).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i too wondered about x::DenseArray. but that is what spmv! inputs. is there a use case for it there? smpv! also does not check the stride of AP. happy to change it here, but would be nice to keep everything in sync.

if 2*length(AP) < N*(N + 1)
throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
end
return spr!(uplo, N, convert(T, α), x, stride(x, 1), AP)
Copy link
Contributor

@sostock sostock Nov 1, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will use the wrong pointer for x if stride(x,1) is negative. However, no BLAS wrappers except for gemv! (edit: and dot) currently handle negative strides correctly, so it may be ok to leave this for a separate PR.

(Support for negative strides in gemv! was added in #41513 and I plan to add it for the other BLAS functions once there is a decision on whether to merge #42012 or #42054.)

stdlib/LinearAlgebra/src/blas.jl Outdated Show resolved Hide resolved
Co-authored-by: Sebastian Stock <42280794+sostock@users.noreply.github.com>
@dkarrasch
Copy link
Member

I think this PR is complete and consistent with the current state of the surrounding functions, so I'll merge once tests pass. I believe we have an issue tracking negative strides, another issue raised was to check the stride for the vector inputs and to catch non-sense uplo arguments. These should be ideally dealt with in a consistent fashion in other PRs.

@bjarthur
Copy link
Contributor Author

@iolaszlo in your PR for spmv you included DenseArray in the method signature. what is the use case for that? for my PR for spr i used your code as a template, and folks are wondering why not just use AbstractVector alone. thx.

@dkarrasch
Copy link
Member

dkarrasch commented Nov 17, 2021

@iolaszlo in your PR for spmv you included DenseArray in the method signature. what is the use case for that? for my PR for spr i used your code as a template, and folks are wondering why not just use AbstractVector alone. thx.

Some discussion can be found here: #34320 (comment). I think there are several places in the BLAS context where we consider any contiguous array simply as a list, no matter what its dimensionality is "officially".

@dkarrasch
Copy link
Member

Some more thoughts. Regarding DenseArrays, that doesn't seem to hurt, and generalizes the methods to any array that is stored contiguously, which is pretty much what BLAS methods need, IIUC. We could test that behavior by calling the new BLAS functions (i) on reshaped AL arrays, and (ii) on AL vectors that have appended, say, ones, because the BLAS function appears to just ignore padded elements at the end (at least, that's what the length check suggests).

@dkarrasch dkarrasch merged commit ae336ab into JuliaLang:master Nov 22, 2021
@bjarthur bjarthur deleted the bja/spr branch November 22, 2021 13:25
@fkastner
Copy link
Contributor

The docstring seems wrong. It says

Update matrix A as αAx*x',[...]

while it should be A + α*x*x', iiuc.

@bjarthur
Copy link
Contributor Author

dang, sorry about that, my bad. please submit a PR to fix it!

@fkastner
Copy link
Contributor

alright, I opened #44303. Hopefully still in time for v1.8

LilithHafner pushed a commit to LilithHafner/julia that referenced this pull request Feb 22, 2022
LilithHafner pushed a commit to LilithHafner/julia that referenced this pull request Mar 8, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants