PARDISO DSS for Block Structured System?

[Precioso__Ali]

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. 

[Alexander_K_Intel2]

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]

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.