Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6977 Discussions

A question on solving many similar matrices with PARDISO

chen146
Beginner
354 Views

Hello,

I have a quick question about solving many structurely similar sparse matrices with PARDISO.

I particular, if I have many linear systems A_i x_i = b_i , where i=1,...,N and A_i have the same sparse structure.

I know that only one pointer array pt is needed for solving all of the systems, and I can solve them by the following pseudo code:

-----------------------------------------------------

do i=1,N

PARDISO phase 1 (A_i ,pt)

PARDISO phase 2 (A_i, pt)

PARDISO phase 3 (A_i,b_i,pt)

end do

-----------------------------------------------------

My question is that if it is OK to pull out (and if it is benefitial to do so) the phase 1 and/or phase 2 of PARDISO out of the do loop. That is,

-----------------------------------------------------

PARDISO phase 1 (A_1 ,pt)

do i=1,N

PARDISO phase 2 (A_i, pt)

PARDISO phase 3 (A_i,b_i,pt)

end do

-----------------------------------------------------

or,

-----------------------------------------------------

PARDISO phase 1 (A_1 ,pt)

PARDISO phase 2 (A_1, pt)

do i=1,N

PARDISO phase 3 (A_i,b_i,pt)

end do

-----------------------------------------------------

I have tried all the three possibility above, and it seems like all of them give the same (correct) result. The computational time needed are also almost the same. This experimental result seems to be weired to me, because my intuition is that only the phase 1 (symbolic factorization) can be done once for all the matrices. The phase 2 (numerical factorization) should be done independently for each matrix. I can not find an explanation nor example in the MKL package/manual to tell me the proper way to do this.

Thank you.

0 Kudos
4 Replies
Alexander_K_Intel2
354 Views

Hi,

Could you provide additional information, like mtype and iparm parameters? I asked because, for example, it is correct to call reorder step once only for several types of matrix or factorization step could be switch on iterative algorithm with precondition when some iparm parameters switch on. Moreover, there is a testcase when reordering step spent more time than factorization and solving step. So we need to have more detailes to investigate your issue.

With best regards,

Alexander Kalinkin

0 Kudos
chen146
Beginner
354 Views

Hello Alexander,

Thank you for the reply. I am actually new to PARDISO, so I am using the iparm and mtype parameters the same as the example file "pardiso_unsym_complex_f.f" in the MKL example sources for my test. That is, mtype=13 (complex, unsym) and iparm is:

      iparm(1) = 1
      iparm(2) = 2
      iparm(3) = 1
      iparm(4) = 0
      iparm(5) = 0
      iparm(6) = 0
      iparm(7) = 0
      iparm(8) = 2
      iparm(9) = 0
      iparm(10) = 13
      iparm(11) = 1
      iparm(12) = 0
      iparm(13) = 1
      iparm(14) = 0
      iparm(15) = 0
      iparm(16) = 0
      iparm(17) = 0
      iparm(18) = -1
      iparm(19) = -1
      iparm(20) = 0

In the real program I am working on, I am trying to solve a complex symmetric shift-invert problem (s*I-A)x=b, where s is a real number, I is the identity matrix, A is a complex symmetric matrix, x and b are vectors. I need to solve this for many different scalers s. Would you mind to give me some suggenstions about how to solve this type of problem efficiently with PARDISO, like what is the good choice of iparm? Thank you and I appreciate it.

Best regards,

Chih-Chieh

0 Kudos
Alexander_K_Intel2
354 Views

Hi,

From my point of view for your case better to use pardiso in following style:

do i=1,N

PARDISO phase 13 (A_i ,pt)

PARDISO phase -1 (pt)

end do

It's possible to use LU decomposition from A_{i} to calculate A_{i+1}, A_{i+2} and etc. But quality of this precondition depend on difference between A_{i} and A_{i+1} so you need to check it. So, it's better to use simple call to pardiso interface with your iparm

With best regards,

Alexander Kalinkin 

0 Kudos
chen146
Beginner
354 Views

Hello Alexander,

Thank you very much. Before I try PARDISO, I have tried to use the LU decomposition of A[0] or A as preconditioners with an iterative solver GMRES (not the one in MKL) and expeiremnts show that the use of preconditioners gives no siginificant improvement in the efficiency. So I think I will go with your suggestion. Thanks again.

Best regards,

Chih-Chieh

0 Kudos
Reply