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.

Sign up for GitHub

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 p-norm functionality. #688

Closed
wants to merge 7 commits into from
Closed

Add p-norm functionality. #688

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9180f52
Add p-norm functionality.
SomTambe Jan 29, 2021
e424e01
Add tests.
SomTambe Feb 1, 2021
1e34822
Set CUDA.allowscalar to false before testing.
SomTambe Feb 1, 2021
fcc1c26
Modified scalar iteration to a macro.
SomTambe Feb 1, 2021
11dd2ac
Changed to if/ else if/ else.
SomTambe Feb 2, 2021
f23d649
Corrected syntax.
SomTambe Feb 2, 2021
c04d918
Corrected syntax again.
SomTambe Feb 2, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions lib/cublas/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,18 @@ end
LinearAlgebra.norm(x::DenseCuArray{<:CublasFloat}) = nrm2(x)
LinearAlgebra.BLAS.asum(x::StridedCuArray{<:CublasFloat}) = asum(length(x), x)

function LinearAlgebra.norm(x::DenseCuArray{<:CublasFloat}, p::Integer)
SomTambe marked this conversation as resolved.
Show resolved Hide resolved
if p==1
return CUBLAS.asum(length(x),x)
end
if p==2
return LinearAlgebra.norm(x)
end
if p>2
return LinearAlgebra.tr(LinearAlgebra.Diagonal(abs.(x))^p)^(1/p)
Copy link
Member

Choose a reason for hiding this comment

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

The problem is that these operations are implemented using CUBLAS, so are themselves restricted to the types CUBLAS supports (plain, basic C types). Ideally we'd have something more generic. What about #84 (comment), is that not valid or slower?

Copy link
Member

Choose a reason for hiding this comment

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

Such a generic versions could also go into GPUArrays.jl, while a specialized version for CUBLAS-supported types could live in CUDA.jl.

Copy link
Author

Choose a reason for hiding this comment

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

@maleadt When we use the method you mentioned in #84, it throws out a scalar indexing warning, which wouldn't work when CUDA.allowscalar(false) is asserted. That is why I have implemented it using CUBLAS based operations.

Copy link
Author

Choose a reason for hiding this comment

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

@maleadt Could you elaborate on the types CUBLAS does not support?

Copy link
Member

Choose a reason for hiding this comment

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

mapreduce is supported by CUDA.jl, I don't understand how that would trigger scalar indexing?

If you look at the CUBLAS docs, you'll see it's a C library that only supports a limited set of element types. That's why we have type unions like CUBLASFloat in CUDA.jl

end
end

function LinearAlgebra.axpy!(alpha::Number, x::StridedCuArray{T}, y::StridedCuArray{T}) where T<:CublasFloat
length(x)==length(y) || throw(DimensionMismatch("axpy arguments have lengths $(length(x)) and $(length(y))"))
axpy!(length(x), alpha, x, y)
Expand Down
9 changes: 9 additions & 0 deletions test/cublas.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
using CUDA.CUBLAS
using CUDA.CUBLAS: band, bandex
using GPUArrays

using LinearAlgebra

Expand Down Expand Up @@ -84,6 +85,14 @@ end
mul!(y, f(A), x, Ts(1), Ts(2))
@test Array(dy) ≈ y
end

@testset "norm" begin
x, y, z = CUDA.rand(elty, 1), CUDA.rand(elty, 2), CUDA.rand(elty, m)
@disallowscalar @test typeof(norm(x,1)) <: Real
@disallowscalar @test typeof(norm(y,2)) <: Real
@disallowscalar @test typeof(norm(z,m)) <: Real
end

@testset "banded methods" begin
# bands
ku = 2
Expand Down