JuliaGPU/CUDA.jl

inv or logdet for symmetric matrices

Open

#977 geöffnet am 14. Juni 2021

Auf GitHub ansehen
 (2 Kommentare) (2 Reaktionen) (0 zugewiesene Personen)Julia (274 Forks)batch import
cuda arrayenhancementgood first issue

Repository-Metriken

Stars
 (1.408 Stars)
PR-Merge-Metriken
 (Durchschn. Merge 5T 5h) (16 gemergte PRs in 30 T)

Beschreibung

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