Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- pardiso forward substitution problems

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted
##

may_ka

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-15-2018
08:11 PM

11 Views

pardiso forward substitution problems

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

3 Replies

Highlighted
##

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-16-2018
03:30 AM

11 Views

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.

Highlighted
##

may_ka

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-16-2018
04:48 AM

11 Views

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.

Highlighted
##

may_ka

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-16-2018
04:58 AM

11 Views

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

For more complete information about compiler optimizations, see our Optimization Notice.