JuliaGPU/CUDA.jl

inv or logdet for symmetric matrices

Open

#977 opened on Jun 14, 2021

View on GitHub
 (2 comments) (2 reactions) (0 assignees)Julia (274 forks)batch import
cuda arrayenhancementgood first issue

Repository metrics

Stars
 (1,408 stars)
PR merge metrics
 (Avg merge 5d 5h) (16 merged PRs in 30d)

Description

If you are looking like me to compute logdet on Symmetric matrices you need to go call some "low-level" functions in CUDA.CUSOLVER This is connected to #116

Here is an unsafe gist for logdet :

function LinearAlgebra.logdet(A::CUDA.CuMatrix)
    d_A = copy(A) # Create a copy
    _, info = CUDA.CUSOLVER.potrfBatched!('L', [d_A]) # Update d_A in-place to be passed to a Cholesky constructor
    L = LinearAlgebra.Cholesky(d_A, 'L', first(info)).L # Create the corresponding lower-triangular matrix
    return 2 * sum(log, diag(L)) # Compute logdet by summing the elements of the diagonal
end

I suppose this could be improved and be more faithful to the implementation in LinearAlgebra but that's further work For inv one can do a similar thing :

function LinearAlgebra.inv(A::CUDA.CuMatrix)
    d_A = copy(A)
    _, info = CUDA.CUSOLVER.potrfBatched!('L', [d_A])
    L = LinearAlgebra.Cholesky(d_A, 'L', first(info)).L
    return (I / L) / L'
end

If you'd rather use lu instead replace the second and third line of each algorithm by

    d_A, p_inv, info = CUDA.CUSOLVER.getrf!(d_A) # Update d_A in-place to be passed to be passed to a LU constructor
    fact = LinearAlgebra.LU(d_A, Int64.(collect(p_inv), zero()) # Create a LU object

and change the last line accordingly

Contributor guide