How does selected inversion work?

[timweiland]

Hi,

Thanks for this great package! I was wondering if and how Pardiso's selected inversion works.
I've followed the process in Pardiso's manual: Call Pardiso in Phase 12, then again in Phase -22 (selected inversion). And it runs through, but where's the actual output (i.e. the sparse inverse)? The manual says that for IPARM(36) = 0 it overwrites the internal L and U factors - can we access these somehow?

I'd be grateful for a small example or any kind of pointers.

[mipals]

Hey there,

I have never done this so I was also a little puzzled as to how this exactly works. From what I understand is that you have to give in a matrix with non-zeros on the entries of the inverse that you aim to compute. The results will then be overwritten in this matrix. I have translated the example pardiso_sym.c below

using SparseArrays
using Pardiso
using LinearAlgebra

ia = [0, 4, 7, 9, 11, 14, 16, 17, 18 ] .+ 1 # C-to-Julia
ja = [0,    2,       5, 6,
                     1, 2,    4,
                        2,             7,
                           3,       6,
                              4, 5, 6,
                                 5,    7,
                                    6,
                                       7] .+ 1 # C-to-Julia
a = [7.0,      1.0,           2.0, 7.0,
                      -4.0, 8.0,           2.0,
                            1.0,                     5.0,
                                 7.0,           9.0,
                                      5.0, -1.0, 5.0,
                                           0.0,      5.0,
                                               11.0,
                                                     5.0];

A = SparseMatrixCSC(8,8,ia,ja,a);

ps = Pardiso.PardisoSolver()
Pardiso.set_matrixtype!(ps,Pardiso.REAL_SYM_INDEF)
Pardiso.ccall_pardisoinit(ps)
ps.msglvl = Pardiso.MESSAGE_LEVEL_ON
b = convert.(Float64,Vector(1:8))
x = similar(b)

# Solve
Pardiso.set_phase!(ps,Pardiso.Phase(13))
ps.iparm[8] = 1
Pardiso.pardiso(ps,x,A,b)
Symmetric(A,:L)*x - b

set_phase!(ps, Pardiso.SELECTED_INVERSION)
ps.iparm[36] = 1
B = similar(A) # This should be the desired sparsity pattern of the inverse
Pardiso.pardiso(ps,x,B,b)

[timweiland]

Thanks for the very fast reply and the useful example! I understand how it works now.
I think it would be super helpful to have a small example in the README, perhaps particularly the last four lines of your example code (that's where I messed up). Let me know if I should make a PR for this.
In any case, thanks a lot!