topic The wrong result when using Pardiso solver to solve linear equations in parallel (OpenMP) in Intel® oneAPI Math Kernel Library

The wrong result when using Pardiso solver to solve linear equations in parallel (OpenMP)

Ewigstern — Tue, 13 Sep 2022 01:41:28 GMT

Hello,

I recently found a strange phenomenon when calling the pardiso solver to solve a linear equation, the details are as follows:

For a linear equation with multiple right-hand sides: AX=B, where A is an n*n matrix, B is an n*m matrix, and m is the number of right-hand sides.

When I solve the equation in the following way (pseudo code), the result is correct.

do i=1,m
call pardiso(A,B(i),X(i))
end do

correct result:

However, when I solve the equation in the following way (pseudo code), all or some of the threads compute the same result, i.e. all or some of the columns of X are the same.

Notice:

1. I have confirmed that the B(i) entered in each loop is correct.
2. These same calculation results may be the correct result of a certain thread, or may be the initial 0 value.
3. In rare cases, all threads can get different calculation results, and the results are correct.

!$OMP Parallel Do Num_Threads(p)
do i=1,m
call pardiso(A,B(i),X(i))
end do
!$OMP end Parallel Do

wrong or unexpected result:

Note: Parallel results are not repeatable (i.e., every time you run the code, you get a different result, and most of the time the results are unexpected). In addtion, parallel results are often wrong or unexpected whether I use sequential or parallel MKL library.

Such a phenomenon is very strange, and it has puzzled me for a long time. I would like to know what is causing this? Is it not possible to call the pardiso solver in the OpenMP way? Or do I need to adjust the input parameters of the pardiso solver when calling the pardiso solver in OpenMP? Any help would be greatly appreciated!

Some information that may be helpful in solving the problem is as follows：

1. MKL Version: Intel(R) oneAPI Math Kernel Library Version 2022.1-Product Build 20220311 for Intel(R) 64 architecture applications.

2. Attachment is the simple test code I used. Note: This test code works in VS Enterprise 2019.

3. OS environment details: Windows 10 v1909, i7 7700K, Visual Studio Enterprise 2019.

PS: Vidya, thank you so much for getting back to me, but for some reason I don't know, I can't get back to you under the original post. I'm so sorry about that, and the only thing I can do is post it again. I hope you can see it. @VidyalathaB_Intel

Thanks,
Ewigstern

Re: The wrong result when using Pardiso solver to solve linear equations in parallel (OpenMP)

VidyalathaB_Intel — Wed, 14 Sep 2022 13:03:22 GMT

Hi Ewigstern,

Thanks for reaching out to us.

Could you please let us know the reason/use case behind calling pardiso in OpenMP parallel region? Because according to the MKL manual

"By default, Intel® oneAPI Math Kernel Library uses the number of OpenMP threads equal to the number of physical cores on the system."

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-windows-developer-guide/top/managing-performance-and-memory/improving-performance-with-threading.html

To do that you can run the code by setting MKL to the /Qmkl: Parallel option.

>>I need to adjust the input parameters of the pardiso solver when calling the pardiso solver in OpenMP

Maybe you can refer to the below link which might help you https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html

Regards,

Vidya.

Re:The wrong result when using Pardiso solver to solve linear equations in parallel (OpenMP)

VidyalathaB_Intel — Wed, 21 Sep 2022 09:10:03 GMT

Hi,

Reminder:

As we haven't back from you, could you please provide us with an update regarding the issue ?

Regards,

Vidya.

Re: The wrong result when using Pardiso solver to solve linear equations in parallel (OpenMP)

Fengrui — Wed, 26 Feb 2025 23:38:46 GMT

Calling oneMKL Pardiso in a parallel region with the same handle pt is not valid usage. Because pt is the handle to some internal data structures that are shared by all the threads. We will add a note about it in the 2025.1 developer reference.