Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
6427 Discussions

PARDISO DSS for Block Structured System?

Precioso__Ali
Beginner
138 Views

I am trying to figure out if PARDISO DSS Interface can be used at the end to solve the following block structured system:

[A11, A12; A21, A22]*[x1,x2]=[b1,b2]

where A11, A12, A21 and A22 are very large sparse matrices and the following Block Gauss-Seidel and Block-Jacobi solutions are needed at the same time :

BJ:       [A11, 0; 0, A22]*[y1,y2]=[b1,b2]

BGS:   [A11, 0; A21, A22]*[y1,y2]=[b1,b2]

For Block Jacobi, DSS Interface can be used to factorize the sparse representation of the diagonal system [A11, 0; 0, A22], to be used repeatedly afterward. However, how can I use the factorization data of the diagonal system [A11, 0; 0, A22] for my Block Gauss-Seidel solution? Unfortunately, the generated  DSS handle is like a black-box, not allowing me to extract computed components from the routine. I need the factorization results to apply the backward-forward solutions partially, once at a time for the rows corresponding to each block, and then update the right-hand side vector accordingly for the BGS operation. Why there is no method to get the L component out of PARDISO routine, same as an already existing pardiso_getdiag for the diagonal component? Any insight on this is appreciated. 

0 Kudos
2 Replies
Alexander_K_Intel2
138 Views

Hi,

L factors inside pardiso store in specific format and convert it in popular csr will significantly decrease performance. However i didn't catch you point - you want to apply same factorization for block Jacobi and Block Gauss-Zeidel, am i right? But in this case their 22 blocks are different, is not it?

Thanks,

Alex

Precioso__Ali
Beginner
138 Views

Hi Alex,

Thank you for your reply. Yes, I would like to pay the cost of matrix factorization only once, and use it for both Block Jacobi and Block Gauss-Seidel operations afterward. This is indeed possible for any block-structured systems in general. For above simple example, BGS should be performed by the following steps:

a)- solve block 1 for A11*y1=b1 --->> y1=inv(A11)*b1

b)- update the right-hand side vector for the next block of the system:  b2up=b2-A21*y1

c)- solve block 2 for A22*y2=b2up  ----> y2=inv(A22)* b2up

In a block Jacobi iteration, you would simply skip step b.

This means for both operations, it is sufficient and enough if we factorize only the single system AD=[A11, 0; 0, A22] representing the diagonal blocks of the original system. By using the L and D info of AD, step a and step c can be done simply by appropriate indexing to an available content like this:

a) - Get access to the content of L and D from row=1 to  row=(size(A11)), then apply backward-forward substitution to get y1=inv(A11)*b1

c) - Get access to the content of L and D from row=size(A11)+1 to row=size(A11)+size(A22), then apply backward-forward substitution to get y2=inv(A22)*b2up

The black-box nature of the PARDISO DSS does not allow to perform step c, because unfortunately there is no support of getting access to L or at least a control to its manipulation at the backward-forward solve step. Is it going to be a solution to this? Or is it possible to command the backward-forward substitution in DSS handle to some specific rows of the matrix?
 
Thanks for your comment.
 
 
 
 

 

Reply