I am introducing the use of MKL PARDISO for the solution of linear systems in some Finite Element legacy code. The 'old' (i.e. 'MKL-less') version of the code writes on a file the values of the determinant and diagonal elements of the (factorized) coefficient matrix. You could reasonably argue that this is not necessarily useful, or that it is even detrimental for very large matrices. However, for historical reasons, I do not have the choice at the moment and the 'MKL version' of the code must provide the same information.
Now that we're done with the context, let us make some considerations:
1) OneAPI MKL PARDISO interface can provide the diagonal (through getdiag()) but not the determinant of the coefficients matrix
2) The DSS interface for PARDISO can provide the determinant but not the diagonal of the factorized matrix
3) Of course, I do not want to factorize the same matrix twice (once using oneAPI-MKL interface and once using DSS interface) to obtain the two pieces of information.
4) For the moment the only solution I came up with is the following (MTYPE=-2, in my case):
- Let A be the coefficient matrix, its factorization is (or should be, if I interpret correctly what PARDISO does) A=PLDL^T (P being a permutation matrix).
- PARDISO returns d=diag(D). So 'almost' OK for the diagonal (please see question c) below).
- We should also be able to retrieve the determinant as det(A)=sign*det(L)*det(D)*det(L^T), where det(L)=det(L^T)=1, det(D) = Product of all the d(i), and sign=+1 if even number of permutations and sign=-1 if odd number of permutations.
Here are my questions:
a) Are considerations 1) - 4) correct?
b) Is there a way to retrieve the number of permutations occurred during factorization phase?
c) the documentation states that the order of the elements returned by getdiag() may not be coherent with the permutation vector (perm) computed in phase 1. Can I conclude that the order of the diagonal elements coincide with the original order of the equations (before factorization, hence before permutation)? Or this is not the case either?
I hope I've made myself sufficiently clear.
Thank you very much for your help!