JuliaGPU/CUDA.jl

inv or logdet for symmetric matrices

Open

#977 aperta il 14 giu 2021

Vedi su GitHub
 (2 commenti) (2 reazioni) (0 assegnatari)Julia (274 fork)batch import
cuda arrayenhancementgood first issue

Metriche repository

Star
 (1408 star)
Metriche merge PR
 (Merge medio 5g 5h) (16 PR mergiate in 30 g)

Descrizione

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

Guida contributor