Skip to content

oneMKL: return LAPACK info from potrf!/getrf!/sytrf! through LinearAlgebra.LAPACK #754

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

Open
michel2323 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

oneMKL: return LAPACK info from potrf!/getrf!/sytrf! through LinearAlgebra.LAPACK #754

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
66 changes: 45 additions & 21 deletions lib/mkl/wrappers_lapack.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,22 +1,34 @@
# LAPACK status of the support library: `computation_error`'s info (> 0, e.g. the
# order of the leading minor that is not positive definite), 0, or -1 when a SYCL
# exception was caught. A SYCL failure is thrown; the factorization info is returned
# by the `LinearAlgebra.LAPACK` methods below, which `cholesky!`/`lu!` destructure.
function chklapackinfo(info, name)
info < 0 && throw(ArgumentError("oneMKL $name failed with a SYCL exception (status $info)"))
return Int64(info)
end

# potrf
for (bname, fname, elty) in ((:onemklSpotrf_scratchpad_size, :onemklSpotrf, :Float32),
(:onemklDpotrf_scratchpad_size, :onemklDpotrf, :Float64),
(:onemklCpotrf_scratchpad_size, :onemklCpotrf, :ComplexF32),
(:onemklZpotrf_scratchpad_size, :onemklZpotrf, :ComplexF64))
@eval begin
function potrf!(uplo::Char,
A::oneStridedMatrix{$elty})
function _potrf!(uplo::Char,
A::oneStridedMatrix{$elty})
chkuplo(uplo)
n = checksquare(A)
lda = max(1, stride(A, 2))

queue = global_queue(context(A), device())
scratchpad_size = $bname(sycl_queue(queue), uplo, n, lda)
scratchpad = oneVector{$elty}(undef, scratchpad_size)
$fname(sycl_queue(queue), uplo, n, A, lda, scratchpad, scratchpad_size)
info = $fname(sycl_queue(queue), uplo, n, A, lda, scratchpad, scratchpad_size)
unsafe_free!(scratchpad)

return A
return A, chklapackinfo(info, "potrf")
end

potrf!(uplo::Char, A::oneStridedMatrix{$elty}) = first(_potrf!(uplo, A))
end
end

Expand Down Expand Up @@ -76,26 +88,31 @@ for (bname, fname, elty) in ((:onemklSsytrf_scratchpad_size, :onemklSsytrf, :Flo
(:onemklCsytrf_scratchpad_size, :onemklCsytrf, :ComplexF32),
(:onemklZsytrf_scratchpad_size, :onemklZsytrf, :ComplexF64))
@eval begin
function sytrf!(uplo::Char,
A::oneStridedMatrix{$elty},
ipiv::oneStridedVector{Int64})
function _sytrf!(uplo::Char,
A::oneStridedMatrix{$elty},
ipiv::oneStridedVector{Int64})
chkuplo(uplo)
n = checksquare(A)
lda = max(1, stride(A, 2))

queue = global_queue(context(A), device())
scratchpad_size = $bname(sycl_queue(queue), uplo, n, lda)
scratchpad = oneVector{$elty}(undef, scratchpad_size)
$fname(sycl_queue(queue), uplo, n, A, lda, ipiv, scratchpad, scratchpad_size)
info = $fname(sycl_queue(queue), uplo, n, A, lda, ipiv, scratchpad, scratchpad_size)
unsafe_free!(scratchpad)

return A, ipiv
return A, ipiv, chklapackinfo(info, "sytrf")
end

function sytrf!(uplo::Char, A::oneStridedMatrix{$elty})
function _sytrf!(uplo::Char, A::oneStridedMatrix{$elty})
n = checksquare(A)
ipiv = oneVector{Int64}(undef, n)
sytrf!(uplo, A, ipiv)
_sytrf!(uplo, A, ipiv)
end

sytrf!(uplo::Char, A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}) =
_sytrf!(uplo, A, ipiv)[1:2]
sytrf!(uplo::Char, A::oneStridedMatrix{$elty}) = _sytrf!(uplo, A)[1:2]
end
end

Expand All @@ -105,23 +122,27 @@ for (bname, fname, elty) in ((:onemklSgetrf_scratchpad_size, :onemklSgetrf, :Flo
(:onemklCgetrf_scratchpad_size, :onemklCgetrf, :ComplexF32),
(:onemklZgetrf_scratchpad_size, :onemklZgetrf, :ComplexF64))
@eval begin
function getrf!(A::oneStridedMatrix{$elty})
function _getrf!(A::oneStridedMatrix{$elty})
m, n = size(A)
ipiv = oneVector{Int64}(undef, min(m, n))
getrf!(A, ipiv)
_getrf!(A, ipiv)
end

function getrf!(A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64})
function _getrf!(A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64})
m,n = size(A)
lda = max(1, stride(A, 2))

queue = global_queue(context(A), device())
scratchpad_size = $bname(sycl_queue(queue), m, n, lda)
scratchpad = oneVector{$elty}(undef, scratchpad_size)
$fname(sycl_queue(queue), m, n, A, lda, ipiv, scratchpad, scratchpad_size)
info = $fname(sycl_queue(queue), m, n, A, lda, ipiv, scratchpad, scratchpad_size)
unsafe_free!(scratchpad)

return A, ipiv
return A, ipiv, chklapackinfo(info, "getrf")
end

getrf!(A::oneStridedMatrix{$elty}) = _getrf!(A)[1:2]
getrf!(A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}) = _getrf!(A, ipiv)[1:2]
end
end

Expand Down Expand Up @@ -608,14 +629,17 @@ end
# LAPACK
for elty in (:Float32, :Float64, :ComplexF32, :ComplexF64)
@eval begin
LinearAlgebra.LAPACK.potrf!(uplo::Char, A::oneStridedMatrix{$elty}) = oneMKL.potrf!(uplo, A)
# the stdlib's return shapes, with the info that `cholesky!`, `lu!` and
# `bunchkaufman!` destructure; `check` (finite input) is not applied on the device
LinearAlgebra.LAPACK.potrf!(uplo::Char, A::oneStridedMatrix{$elty}) = oneMKL._potrf!(uplo, A)
LinearAlgebra.LAPACK.potrs!(uplo::Char, A::oneStridedMatrix{$elty}, B::oneStridedVecOrMat{$elty}) = oneMKL.potrs!(uplo, A, B)
LinearAlgebra.LAPACK.sytrf!(uplo::Char, A::oneStridedMatrix{$elty}) = oneMKL.sytrf!(uplo, A)
LinearAlgebra.LAPACK.sytrf!(uplo::Char, A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}) = oneMKL.sytrf!(uplo, A, ipiv)
LinearAlgebra.LAPACK.sytrf!(uplo::Char, A::oneStridedMatrix{$elty}) = oneMKL._sytrf!(uplo, A)
LinearAlgebra.LAPACK.sytrf!(uplo::Char, A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}) = oneMKL._sytrf!(uplo, A, ipiv)
LinearAlgebra.LAPACK.geqrf!(A::oneStridedMatrix{$elty}) = oneMKL.geqrf!(A)
LinearAlgebra.LAPACK.geqrf!(A::oneStridedMatrix{$elty}, tau::oneStridedVector{$elty}) = oneMKL.geqrf!(A, tau)
LinearAlgebra.LAPACK.getrf!(A::oneStridedMatrix{$elty}) = oneMKL.getrf!(A)
LinearAlgebra.LAPACK.getrf!(A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}) = oneMKL.getrf!(A, ipiv)
LinearAlgebra.LAPACK.getrf!(A::oneStridedMatrix{$elty}; check::Bool = true) = oneMKL._getrf!(A)
LinearAlgebra.LAPACK.getrf!(A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}; check::Bool = true) =
oneMKL._getrf!(A, ipiv)
LinearAlgebra.LAPACK.getrs!(trans::Char, A::oneStridedMatrix{$elty}, ipiv::oneStridedVector{Int64}, B::oneStridedVecOrMat{$elty}) = oneMKL.getrs!(trans, A, ipiv, B)
LinearAlgebra.LAPACK.ormqr!(side::Char, trans::Char, A::oneStridedMatrix{$elty}, tau::oneStridedVector{$elty}, C::oneStridedVecOrMat{$elty}) = oneMKL.ormqr!(side, trans, A, tau, C)
LinearAlgebra.LAPACK.orgqr!(A::oneStridedMatrix{$elty}, tau::oneStridedVector{$elty}) = oneMKL.orgqr!(A, tau)
Expand Down
50 changes: 50 additions & 0 deletions test/onemkl.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1893,6 +1893,38 @@ end
@test B ≈ collect(d_B)
end

@testset "LAPACK.potrf! info -- cholesky" begin
A = rand(elty, n, n)
A = A*A' + n*I
d_A, info = LAPACK.potrf!('L', oneArray(A))
@test info == 0
@test tril(collect(d_A)) ≈ tril(LAPACK.potrf!('L', copy(A))[1])
C = cholesky(Hermitian(oneArray(A), :L))
@test issuccess(C)
@test tril(collect(C.factors)) ≈ cholesky(Hermitian(A, :L)).L

# not positive definite at column j: info == j, as LAPACK
N = 100
G = rand(elty, N, N)
H = Matrix(Hermitian(G*G' + N*I))
for j in (1, 17, 32)
Hj = copy(H); Hj[j, j] = -one(elty)
@test LAPACK.potrf!('L', oneArray(Hj))[2] == j
@test LAPACK.potrf!('U', oneArray(Hj))[2] == j
@test !issuccess(cholesky(Hermitian(oneArray(Hj), :L); check = false))
end
# oneMKL (2026.1, Intel Data Center GPU Max) names the failing column relative to
# its 32-wide panel beyond column 32 (column 50 comes back as 18, or as 1042 for
# complex types); the factor's diagonal still holds the failed pivot at column j
for j in (33, 50, N)
Hj = copy(H); Hj[j, j] = -one(elty)
d_Hj, info = LAPACK.potrf!('L', oneArray(Hj))
@test_broken info == j
d = real.(diag(collect(d_Hj)))
@test findfirst(x -> !(isfinite(x) && x > 0), d) == j
end
end

# @testset "sytrf!" begin
# A = rand(elty,n,n)
# A = A + A'
Expand All @@ -1917,6 +1949,24 @@ end
@test h_A ≈ Array(d_A)
end

@testset "LAPACK.getrf! info -- lu" begin
A = rand(elty, m, m) + m*I
d_A, d_ipiv, info = LAPACK.getrf!(oneArray(A))
h_A, ipiv, h_info = LAPACK.getrf!(copy(A))
@test info == h_info == 0
@test collect(d_A) ≈ h_A
@test collect(d_ipiv) == ipiv
F = lu(oneArray(A))
@test issuccess(F)
@test collect(F.factors) ≈ h_A
@test collect(F.ipiv) == ipiv

# exactly singular: info is the zero pivot, as LAPACK
S = copy(A); S[:, 3] .= 0
@test LAPACK.getrf!(oneArray(S))[3] == LAPACK.getrf!(copy(S))[3] == 3
@test !issuccess(lu(oneArray(S); check = false))
end

@testset "getrf_batched! -- getri_batched!" begin
bA = [rand(elty, m, m) for i in 1:p]
d_bA = oneMatrix{elty}[]
Expand Down
Loading