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

[CUSOLVER] Add new methods for \ and inv #2117

Merged
merged 4 commits into from
Oct 18, 2023
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
Fix the tests with complex numbers
  • Loading branch information
amontoison committed Oct 18, 2023
commit 7b444cf816f7b7884f536a8d0fcd34af69183beb
44 changes: 17 additions & 27 deletions lib/cusolver/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -65,20 +65,15 @@ function Base.:\(_A::CuMatOrAdj, _B::CuOrAdj)
return X
end

for (wrapper, blas_types) in ((:Symmetric, :BlasReal),
(:Hermitian, :BlasFloat))
@eval begin
function Base.:\(A::$wrapper{<:$blas_types,<:CuMatOrAdj}, _B::CuOrAdj)
uplo = A.uplo
A, B = copy_cublasfloat(_A.data, _B)

# LDLᴴ decomposition with partial pivoting
factors, ipiv, info = sytrf!(uplo, F)
ipiv = CuVector{Int64}(ipiv)
X = sytrs!(uplo, factors, ipiv, B)
return X
end
end
function Base.:\(_A::Symmetric{<:Any,<:CuMatOrAdj}, _B::CuOrAdj)
uplo = A.uplo
A, B = copy_cublasfloat(_A.data, _B)

# LDLᴴ decomposition with partial pivoting
factors, ipiv, info = sytrf!(uplo, A)
ipiv = CuVector{Int64}(ipiv)
X = sytrs!(uplo, factors, ipiv, B)
return X
end

# patch JuliaLang/julia#40899 to create a CuArray
Expand Down Expand Up @@ -370,19 +365,14 @@ function LinearAlgebra.inv(A::StridedCuMatrix{T}) where T <: BlasFloat
return B
end

for (wrapper, blas_types) in ((:Symmetric, :BlasReal),
(:Hermitian, :BlasFloat))
@eval begin
function LinearAlgebra.inv(A::$wrapper{T,<:StridedCuMatrix{T}}) where T <: $blas_types
n = checksquare(A)
F = copy(A.data)
factors, ipiv, info = sytrf!(A.uplo, F)
ipiv = CuVector{Int64}(ipiv)
B = CuMatrix{T}(I, n, n)
sytrs!(A.uplo, factors, ipiv, B)
return B
end
end
function LinearAlgebra.inv(A::Symmetric{T,<:StridedCuMatrix{T}}) where T <: BlasFloat
n = checksquare(A)
F = copy(A.data)
factors, ipiv, info = sytrf!(A.uplo, F)
ipiv = CuVector{Int64}(ipiv)
B = CuMatrix{T}(I, n, n)
sytrs!(A.uplo, factors, ipiv, B)
return B
end

for (triangle, uplo, diag) in ((:LowerTriangular, 'L', 'N'),
Expand Down
39 changes: 15 additions & 24 deletions test/libraries/cusolver/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,14 @@ k = 1
@test Array(dI) ≈ I
end

@testset "inv -- symmetric -- hermitian" begin
@testset "inv -- symmetric" begin
A = rand(elty,n,n)
A = A + A'
dA = Hermitian(CuArray(A))
A = A + tranpose(A)
dA = Symmetic(CuArray(A))
dA⁻¹ = inv(dA)
dI = dA.data * dA⁻¹
display(dI)
@test Array(dI) ≈ I

if elty <: Real
dA = Symmetric(CuArray(A))
dA⁻¹ = inv(dA)
dI = dA.data * dA⁻¹
@test Array(dI) ≈ I
end
end

@testset "inv -- triangular" begin
Expand All @@ -39,7 +33,7 @@ k = 1
A = rand(elty,n,n)
A = uplo == 'L' ? tril(A) : triu(A)
A = diag == 'N' ? A : A - Diagonal(A) + I
dA = CuArray(A)
dA = triangle(CuArray(A))
dA⁻¹ = inv(dA)
dI = dA.data * dA⁻¹
@test Array(dI) ≈ I
Expand Down Expand Up @@ -712,28 +706,25 @@ end
], elty2 in [
Float16, Float32, Float64, ComplexF16, ComplexF32, ComplexF64, Int32, Int64, Complex{Int32}, Complex{Int64}
]
@testset "Square symmetric and hermitian linear systems" begin
@testset "Symmetric linear systems" begin
A = rand(elty1,n,n)
A = A + A'
A = A + transpose(A)
B = rand(elty2,n,5)
b = rand(elty2,n)
d_A = CuArray(A)
d_B = CuArray(B)
d_b = CuArray(b)
cublasfloat = promote_type(Float32, promote_type(elty1, elty2))
for wrapper in (Symmetric, Hermitian)
(wrapper == Symmetric) && !(cublasfloat <: Real) && continue
Af = wrapper(cublasfloat.(A))
Bf = cublasfloat.(B)
bf = cublasfloat.(b)
@test Array(d_A \ d_B) ≈ (Af \ Bf)
@test Array(d_A \ d_b) ≈ (Af \ bf)
@inferred d_A \ d_B
@inferred d_A \ d_b
end
Af = Symmetric(cublasfloat.(A))
Bf = cublasfloat.(B)
bf = cublasfloat.(b)
@test Array(d_A \ d_B) ≈ (Af \ Bf)
@test Array(d_A \ d_b) ≈ (Af \ bf)
@inferred d_A \ d_B
@inferred d_A \ d_b
end

@testset "Square unsymmetric linear systems" begin
@testset "Square and unsymmetric linear systems" begin
A = rand(elty1,n,n)
B = rand(elty2,n,5)
b = rand(elty2,n)
Expand Down
2 changes: 1 addition & 1 deletion test/libraries/cusolver/dense_generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ p = 5
for uplo in ('L', 'U')
A = rand(elty,n,n)
B = rand(elty,n,p)
A = A + A'
A = A + transpose(A)
d_A = CuMatrix(A)
d_B = CuMatrix(B)
d_A, d_ipiv, _ = CUSOLVER.sytrf!(uplo, d_A)
Expand Down