What did I do wrong with LAPACK.potri!?

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.

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.