Showing results for 
Search instead for 
Did you mean: 

Help with the Intel MKL interface to PARDISO

Dear Intel MKL experts, 

I would like to ask for help with the use
of the Intel MKL interface to PARDISO.

Assume that you have a sparse matrix
A partitioned into 2x2 blocks

A= [ A_II A_IG ]
[ A_GI A_GG ],

and that you want to solve linear systems
with A and A_II using only one sparse LU
decomposition of the global system A:

A = [L_II 0] [U_II U_IG]
[L_GI L_GG] [0 U_GG]

Neglecting pivoting for numerical
stability, the solution of linear systems
of the form A_II x_I = b_I can be performed
as follows:

1) Solve [L_II 0] [y_I] = [b_I]
[L_GI L_GG] [y_G] = [ * ]

2) Y_G <- 0

3) Solve [U_II U_IG] [ x_I ] = [y_I]
[0 U_GG] [ x_G=0 ] = [ 0 ]

This algorithm can be implemented with
the Intel MKL interface of PARDISO because
it allows to separate among forward and
backward substitution, so that you can
zero-out Y_G between steps 1 and 3. The problem
that I have is that this workaround
is no longer valid when numerical pivoting enters the scene
(for example, with symmetric indefinite matrices)

Do you know whether there is a simple strategy
to tell PARDISO to restrict numerical pivoting
to a particular block (i.e., A_II) or otherwise
to tell PARDISO to solve a subblock given a
global LU factorization of the global system ?.
I want to avoid to compute both a factorization
of A_II and A.

Thanks a lot for your help.

Best regards,
0 Kudos
1 Reply

Hi Alberto,
Pardiso is direct solver based on LU decomposition with global reordering of initial matrix. In your case you want to set specific reordering to PARDISO so it could be implemented by 2 variants:
1. Use own perm vector by set iparm[4]=1
2. Use partial solve by set iparm[30] = 1
Nevertheless using this parameters could decrease performance and increase internal memory size of PARDISO so I strongly recommend you just factorize by pardiso 2 different matrix: A_II and full initial matrix.
With best regards,
Alexander Kalinkin