pardiso forward substitution problems

[may_ka]

Hi all,

I have difficulties to line up results from pardiso forward substitution with results from explicit trials (eg. in R or pyhton).

Let "C" be a positive definite symmetric squared matrix and "v" be a vector of 1 with length equal to the rows of "C".

When solving "Cb=v" pardiso results looking good with phase 33, where "b" contains the sum over rows of "C^{-1}".

Now "C" maybe decomposed into "LL' " (cholesky factor). A forward substitution would solve "Lb=v", thus the result must be the sum over the rows of "L^{-1}".

However, when feeding the same matrix "C" into pardiso, but setting phase to 331 the results doesn't seem to reflect the above. I guessed that it might be due to permutation. Retrieving the permutation by setting iparm(5) to reshuffle the result vector did not do anything. Moreover, since the range fo the results from pardiso where different to an explicit operation (taking the row sums of "L^{-1}"), the cause cannot the permutation.

Did I miss anything?? For instance is pardiso's LL' not the cholesky of C??

Any advice appreciated.

Thanks

NB: iparm(36) is zero

[mecej4]

Cb = v and Lb = v are in conflict, since the same b cannot be the solution of the full system as well as the triangular system. Had that been possible, there would have been no need for the back-substitution part of the Cholesky algorithm.

"Forward Substitution"? The usual description uses "Forward Elimination" and "Back Substitution"; see, for example, https://en.wikipedia.org/wiki/Gaussian_elimination .

Please rewrite your argument using consistent and unambiguous notation.

[may_ka]

I think I could help myself.

From what I understand is pardiso is doing LL'=PCP', where P is a permutation matrix. I retrieved L from pardiso and observed that it is not a lower triangular matrix. However, the retrieved L yield C when doing LL'. Therefore it appears as if the retrieved L is actually P^{-1}L.

[may_ka]

Hi

thanks for unambiguous comment.

With regard to "substitution" vs "elimination" you may go after the manual authors as well. Just a small read from page 2495 of the 2017 manual:

If it is set to MKL_DSS_FORWARD_SOLVE, the forward substitution (corresponding to phase = 331 in Intel
MKL PARDISO) is performed;
if it is set to MKL_DSS_DIAGONAL_SOLVE, the diagonal substitution (corresponding to phase = 332 in Intel
MKL PARDISO) is performed, if possible;
if it is set to MKL_DSS_BACKWARD_SOLVE, the backward substitution (corresponding to phase = 333 in Intel
MKL PARDISO) is performed.
For more details about using these substitutions for different types of matrices, see Separate Forward and
Backward Substitution in the Intel MKL PARDISO solver description.

Moreover, of course Cb=v and Lb=v are not in conflict. While "v" is constant the content of "b" changes, from the row sum of C^{-1} to the row sum of L^{-1}. This is written explicitly in the initial thread.

I appreciate comments, but I even more appreciate if the commentator has read the question.

Thanks a lot.

Cheers