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 for the diagonal component? Any insight on this is appreciated.
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?
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