Add spmv! to BLAS in stdlib/LinearAlgebra #34320

iolaszlo · 2020-01-08T23:48:36Z

No description provided.

dkarrasch

This looks great. Just one comment on the allowed signature, with relevance for the previous hpmv PR.

stdlib/LinearAlgebra/src/blas.jl

dkarrasch

I think my previous concern on "enforcing 1-dimensionality" was void, so that's resolved. I'd suggest though to explain the packed format explicitly, since the link to the LAPACK page is rather indirect, one has to search on their webpage to get to the specific page. This might carry over to the hpmv code.

stdlib/LinearAlgebra/src/blas.jl

dkarrasch · 2020-01-10T14:34:31Z

stdlib/LinearAlgebra/src/blas.jl

+    if length(AP) < Int64(N*(N+1)/2)
+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    GC.@preserve x y AP spmv!(uplo, N, convert(T, α), AP, pointer(x), BlasInt(stride(x, 1)), convert(T, β), pointer(y), BlasInt(stride(y, 1)))


Do we need to preserve AP explicitly? I guess it's referenced directly and not in danger of being gc'ed.

Maybe I'm too cautious here due to lack of deep GC knowledge, so looking at the other @GC.preserve statements in blas.jl I do not see any crucial difference and reason to drop AP from the preserve statement. Could you maybe explain this a bit more or point out the the thing that I'm missing?

IIUC, then only those are listed which are wrapped by a pointer(...). You could double-check if that's the case, but that was my understanding.

IIUC, then only those are listed which are wrapped by a pointer(...). You could double-check if that's the case, but that was my understanding.

So if I were to drop the pointer(...) around x and y as well then the whole @GC.preserve might be unnecessary? If true this might be tempting to apply to have a more clean API. :)

Possibly following here a pattern similar to that in the symv! implementation strengthens your reading so removing the @GC.preserve alltogether might be good to go. If you are also okay with it I will apply this modification as well.

On the other hand maybe the following is a non-issue but wherever we have a @GC.preserve we also have Union{DenseArray,...} or Union{Unitrange,...} or Union{.., ...}. Maybe these are unrelated but they do appear together so the removal of @GC.preserve would induce to suggest to remove pointer(...) and GC.preserve from all those functions. I don't know whether this is wise as people might have put it this way for a reason unknown to me at the moment.

Probably there is no harm in keeping AP here, so unless someone comes by and helps in resolving this aspect, this should not keep us from merging this.

I am planning other similar pull requests. Would it be okay in your opinion to ask someone within this comment thread whether to use the @GC.preserve to ease further development?

It's better just to pass the arrays, see #34783.

stdlib/LinearAlgebra/test/blas.jl

iolaszlo · 2020-01-10T19:58:33Z

I think my previous concern on "enforcing 1-dimensionality" was void, so that's resolved. I'd suggest though to explain the packed format explicitly, since the link to the LAPACK page is rather indirect, one has to search on their webpage to get to the specific page. This might carry over to the hpmv code.

I agree that part of the docstring changes and more nice tests are necessary and would like to apply them to spmv! and hpmv! as well. Should I apply the spmv! changes in this pull request and start a separate pull requiest for the hpmv! update to keep at least the current pull request atomic in some sense? Or let's try to merge this one and fix both in a separate pull request?

dkarrasch

I'd suggest we work on this until spmv! is perfect, and then edit hpmv! accordingly in another PR. I believe there are only mild adoptions necessary, like the restriction to Real, and some small work on the docstrings, so there is an end in sight. We can raise the GC.gc issue in the next PR again, unless someone comes by and helps out with expertise here.

stdlib/LinearAlgebra/src/blas.jl

dkarrasch · 2020-01-11T11:50:52Z

stdlib/LinearAlgebra/src/blas.jl

+    if length(AP) < Int64(N*(N+1)/2)
+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    GC.@preserve x y AP spmv!(uplo, N, convert(T, α), AP, pointer(x), BlasInt(stride(x, 1)), convert(T, β), pointer(y), BlasInt(stride(y, 1)))


Probably there is no harm in keeping AP here, so unless someone comes by and helps in resolving this aspect, this should not keep us from merging this.

dkarrasch

This is shaping up nicely, just some minor polishing and this may serve as the role model for a minor hpmv! "fix".

dkarrasch · 2020-01-14T08:18:19Z

stdlib/LinearAlgebra/src/blas.jl

+                      (:sspmv_, :Float32))
+    @eval begin
+        function spmv!(uplo::AbstractChar,
+                       n::BlasInt,


I guess we kind of figured out in #34360 that this should be Integer.

dkarrasch · 2020-01-14T08:26:25Z

stdlib/LinearAlgebra/src/blas.jl

+    if length(AP) < Int64(N*(N+1)/2)
+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    GC.@preserve x y AP spmv!(uplo, BlasInt(N), convert(T, α), AP, pointer(x), BlasInt(stride(x, 1)), convert(T, β), pointer(y), BlasInt(stride(y, 1)))


... and here no BlasInt(N) etc. conversions

dkarrasch · 2020-01-14T08:30:37Z

stdlib/LinearAlgebra/test/blas.jl

+                # Define the inputs and outputs of spmv!, y = α*A*x+β*y
+                α = rand(elty)
+                M = rand(elty, n, n)
+                A = Symmetric(M)


I would write

AL = Symmetric(M, :L) AU = Symmetric(M, :U)

...

dkarrasch · 2020-01-14T08:31:04Z

stdlib/LinearAlgebra/test/blas.jl

+                β = rand(elty)
+                y = rand(elty, n)
+
+                y_result_julia = α*A*x+β*y


then

y_result_julia = α*AL*x+β*y

...

dkarrasch · 2020-01-14T08:32:23Z

stdlib/LinearAlgebra/test/blas.jl

+                AP = typeof(A[1,1])[]
+                for j in 1:n
+                    for i in j:n
+                        push!(AP, A[i,j])


...

push!(AP, M[i,j])

and AP = typeof(M[1,1])[] (even though that doesn't matter at all, just for clarity where we take types from)...

dkarrasch · 2020-01-14T08:33:37Z

stdlib/LinearAlgebra/test/blas.jl

+                @test y_result_julia≈y_result_blas_lower
+
+                # Create upper triangular packing of A
+                AP = typeof(A[1,1])[]


... and analogous changes here, with a recomputed

y_result_julia = α*AU*x+β*y

stdlib/LinearAlgebra/src/blas.jl

dkarrasch · 2020-01-15T08:00:14Z

@yuyichao Could you please help us and take a look at the GC.@preserve line. We are unsure whether AP needs to be listed there. It would be great if you could share some considerations regarding GC.@preserve for the future.

@iolaszlo Could you maybe rebase your branch on current master. That might help to resolve some of the (otherwise unrelated) test failures.

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

dkarrasch · 2020-02-17T09:06:54Z

It would be great to get this PR in, and to learn about how to use GC.@preserve correctly for future PRs. Anyone care to advise us? @yuyichao @andreasnoack @ViralBShah ? The remaining question is which variables need to be preserved from gc.

ViralBShah · 2020-02-17T11:30:11Z

+1 to get the PR in, but we need someone else to chime in on GC.@preserve

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

dkarrasch

Thanks for your stamina @iolaszlo! 😃

dkarrasch · 2020-02-19T06:53:44Z

On win32, the LinearAlgebra/blas tests pass, macos64 is a download failure.

andreasnoack

Just a few minor comments, but looks good to me other than that.

stdlib/LinearAlgebra/src/blas.jl

andreasnoack · 2020-02-19T07:20:25Z

stdlib/LinearAlgebra/src/blas.jl

+        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
+    end
+    spmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1))
+    y


...then you can leave it out here. I prefer explicit return spmv!... but it's not required.

andreasnoack · 2020-02-19T07:20:53Z

stdlib/LinearAlgebra/test/blas.jl

+
+                y_result_blas_upper = copy(y)
+                BLAS.spmv!('U', α, AP, x, β, y_result_blas_upper)
+                @test y_result_julia_upper≈y_result_blas_upper


Suggested change

@test y_result_julia_upper≈y_result_blas_upper

@test y_result_julia_upper ≈ y_result_blas_upper

andreasnoack · 2020-02-19T07:21:17Z

stdlib/LinearAlgebra/test/blas.jl

+                @test y_result_julia_lower≈y_result_blas_lower
+
+
+                y_result_julia_upper = α*AU*x+β*y


Suggested change

y_result_julia_upper = α*AU*x+β*y

y_result_julia_upper = α*AU*x + β*y

andreasnoack · 2020-02-19T07:21:27Z

stdlib/LinearAlgebra/test/blas.jl

+                β = rand(elty)
+                y = rand(elty, n)
+
+                y_result_julia_lower = α*AL*x+β*y


Suggested change

y_result_julia_lower = α*AL*x+β*y

y_result_julia_lower = α*AL*x + β*y

andreasnoack · 2020-02-19T07:22:32Z

stdlib/LinearAlgebra/src/blas.jl

+    if N != length(y)
+        throw(DimensionMismatch("x has length $(N), but y has length $(length(y))"))
+    end
+    if length(AP) < Int64(N*(N+1)/2)


Suggested change

if length(AP) < Int64(N*(N+1)/2)

if 2*length(AP) < N*(N + 1)

dkarrasch · 2020-02-21T08:37:49Z

CI failures are in Distributed and Sockets, respectively. I'd say unrelated.

iolaszlo · 2020-02-21T18:56:27Z

@dkarrasch Thank you for the constructive review process.

iolaszlo force-pushed the extend-blas-with-spmv branch from 496bb08 to 42577b2 Compare January 8, 2020 23:51

ViralBShah requested a review from dkarrasch January 9, 2020 00:22

dkarrasch requested changes Jan 9, 2020

View reviewed changes

stdlib/LinearAlgebra/src/blas.jl Outdated Show resolved Hide resolved

dkarrasch requested changes Jan 10, 2020

View reviewed changes

dkarrasch added the linear algebra Linear algebra label Jan 11, 2020

dkarrasch reviewed Jan 11, 2020

View reviewed changes

iolaszlo force-pushed the extend-blas-with-spmv branch from 42577b2 to 116ed29 Compare January 13, 2020 22:58

iolaszlo requested a review from dkarrasch January 13, 2020 23:01

dkarrasch mentioned this pull request Jan 14, 2020

Convert hpmv input to a BlasInt. #34360

Merged

dkarrasch requested changes Jan 14, 2020

View reviewed changes

dkarrasch reviewed Jan 14, 2020

View reviewed changes

stdlib/LinearAlgebra/src/blas.jl Outdated Show resolved Hide resolved

iolaszlo force-pushed the extend-blas-with-spmv branch from 116ed29 to 76dfed6 Compare January 14, 2020 18:45

iolaszlo requested a review from dkarrasch January 14, 2020 18:47

dkarrasch approved these changes Jan 15, 2020

View reviewed changes

iolaszlo force-pushed the extend-blas-with-spmv branch 2 times, most recently from 62044c6 to e6828c8 Compare January 19, 2020 18:21

iolaszlo added a commit to iolaszlo/julia that referenced this pull request Jan 19, 2020

Minor updates to BLAS/hpmv! as in JuliaLang#34320.

375fb39

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

iolaszlo added a commit to iolaszlo/julia that referenced this pull request Jan 19, 2020

Minor updates to BLAS/hpmv! as in JuliaLang#34320. [JuliaLang#34439]

cb11ac7

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

iolaszlo added a commit to iolaszlo/julia that referenced this pull request Jan 19, 2020

Minor updates to BLAS/hpmv! as in JuliaLang#34320. [JuliaLang#34439]

45b2eec

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

andreasnoack mentioned this pull request Feb 17, 2020

Pass arrays instead of pointers when possible in blas.jl #34783

Merged

iolaszlo force-pushed the extend-blas-with-spmv branch from e6828c8 to a602df6 Compare February 18, 2020 22:44

iolaszlo requested review from dkarrasch and andreasnoack February 18, 2020 22:45

iolaszlo added a commit to iolaszlo/julia that referenced this pull request Feb 18, 2020

Minor updates to BLAS/hpmv! as in JuliaLang#34320. [JuliaLang#34439]

bee1d09

Follow the syntax, docstring and testing conventions elaborated in JuliaLang#34320.

dkarrasch approved these changes Feb 19, 2020

View reviewed changes

andreasnoack reviewed Feb 19, 2020

View reviewed changes

Add spmv! to BLAS in stdlib/LinearAlgebra

947cd03

iolaszlo force-pushed the extend-blas-with-spmv branch from a602df6 to 947cd03 Compare February 20, 2020 19:52

iolaszlo requested review from andreasnoack, ViralBShah and dkarrasch February 20, 2020 19:53

andreasnoack approved these changes Feb 20, 2020

View reviewed changes

dkarrasch approved these changes Feb 20, 2020

View reviewed changes

dkarrasch merged commit d8b2209 into JuliaLang:master Feb 21, 2020

birm pushed a commit to birm/julia that referenced this pull request Feb 22, 2020

Add spmv! to BLAS in stdlib/LinearAlgebra (JuliaLang#34320)

a9ed42e

KristofferC pushed a commit that referenced this pull request Apr 11, 2020

Add spmv! to BLAS in stdlib/LinearAlgebra (#34320)

e0f5238

bjarthur mentioned this pull request Nov 16, 2021

add spr! to LinearAlgebra.BLAS #42830

Merged

dkarrasch mentioned this pull request Sep 25, 2023

Implement packed format for symmetric and triangular matrices #51431

Open

	@test y_result_julia_upper≈y_result_blas_upper
	@test y_result_julia_upper ≈ y_result_blas_upper

		@test y_result_julia_lower≈y_result_blas_lower


		y_result_julia_upper = αAUx+β*y

	y_result_julia_upper = αAUx+β*y
	y_result_julia_upper = αAUx + β*y

	y_result_julia_lower = αALx+β*y
	y_result_julia_lower = αALx + β*y

Add spmv! to BLAS in stdlib/LinearAlgebra #34320

Add spmv! to BLAS in stdlib/LinearAlgebra #34320

Conversation

iolaszlo commented Jan 8, 2020

dkarrasch left a comment

Choose a reason for hiding this comment

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iolaszlo Jan 10, 2020 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iolaszlo commented Jan 10, 2020 • edited

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch commented Jan 15, 2020 • edited

dkarrasch commented Feb 17, 2020

ViralBShah commented Feb 17, 2020

dkarrasch left a comment

Choose a reason for hiding this comment

dkarrasch commented Feb 19, 2020

andreasnoack left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch commented Feb 21, 2020

iolaszlo commented Feb 21, 2020

iolaszlo Jan 10, 2020 •

edited

iolaszlo commented Jan 10, 2020 •

edited

dkarrasch commented Jan 15, 2020 •

edited