kai
1
Somehow potri is not giving me the inverse. What did I do wrong? Here is an example:
using LinearAlgebra, Random
Random.seed!(1)
x = rand(5, 5)
V = x'x
LinearAlgebra.copytri!(V, 'U')
Vcopy = rand(5, 5)
copyto!(Vcopy, V)
Then
inv(V) * Vcopy
gives back identity.
But
LinearAlgebra.LAPACK.potri!('U', V)
LinearAlgebra.copytri!(V, 'U')
V * Vcopy
gives something far from identity.
kai
2
Turns out that I need to call potrf
before potri
.
The documentation
Computes the inverse of positive-definite matrix A
after calling potrf!
to find its (upper if uplo = U
, lower if uplo = L
) Cholesky decomposition.
made me think that potri
would call potrf
internally. But that is not the case.
So
LinearAlgebra.LAPACK.potrf!('U', V)
LinearAlgebra.LAPACK.potri!('U', V)
LinearAlgebra.copytri!(V, 'U')
V * Vcopy
would give back identity.
2 Likes
I am not sure why you need to call the LAPACK functions directly. I think that
inv(cholesky(V))
does what you want.