JuliaGPU/CUDA.jl

inv or logdet for symmetric matrices

Open

#977 aberto em 14 de jun. de 2021

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

Métricas do repositório

Stars
 (1.408 stars)
Métricas de merge de PR
 (Mesclagem média 5d 5h) (16 fundiu PRs em 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

Guia do colaborador