topic Re: Pardiso got the wrong Schur complement in Intel® oneAPI Math Kernel Library

Pardiso got the wrong Schur complement

BrianZHANG — Sat, 31 Aug 2024 13:44:06 GMT

Hi, I am calling Pardiso in Fortran to compute the schur complement of sparse matrix. When the matrix is small, I can get the right answer. But I got wrong results when the matrix dimension increases. The mkl I use is the Linux version 2024.2

Here are some of my Pardiso parameters:

Type(MKL_CLUSTER_SPARSE_SOLVER_HANDLE) :: Pt(64) Integer :: Maxfct, Mnum, Mtype, Phase, Neqns, Nrhs, Msglvl, Error Integer :: Perm(Neqns) Integer :: Iparm(64) Do i = 1, 64 Pt(i)%dummy = 0 End Do Maxfct = 1 Mnum = 1 Perm = 0 Nrhs = 1 Error = 0 Msglvl = 0 Mtype = 11 Iparm = 0 Iparm(1) = 1 Iparm(2) = 2 Iparm(3) = OMP_GET_NUM_PROCS() Iparm(8) = 5 Iparm(10) = 13 Iparm(13) = 1 Iparm(21) = 1 Iparm(24) = 1 Iparm(36) = -1 Perm(numLocalInternalFreedomDegrees + 1:) = 1 Phase = 11 Call Pardiso( Pt, Maxfct, Mnum, Mtype, Phase, Neqns, CSR_Matrix_Value, CSR_Row_Index, CSR_Matrix_Column, Perm, Nrhs, Iparm, Msglvl, Right_Vector, Computed_Value, Error ) numNoneZero = Iparm(36) Allocate(KS_Row_Index(numLocalSharedFreedomDegrees + 1)) Allocate(KS_Matrix_Column(numNoneZero), KS_Matrix_Value(numNoneZero)) Call pardiso_export(Pt, KS_Matrix_Value, KS_Row_Index, KS_Matrix_Column, 1, Iparm, Error) Iparm(36) = -1 Phase = 22 Call Pardiso( Pt, Maxfct, Mnum, Mtype, Phase, Neqns, CSR_Matrix_Value, CSR_Row_Index, CSR_Matrix_Column, Perm, Nrhs, Iparm, Msglvl, Right_Vector, Computed_Value, Error )

I have set Iparm(36) to 1 and 2 to calculate dense Schur complement, and the results are not any different from the CSR format results obtained with -1 and -2. I use Matlab to calculate schur complement. In the figure, the blue one is the Matlab's results, and the red one is the Pardiso's results. Pardiso's results are obviously wrong.

Is something wrong in my iparm parameters, or is there a bug in the mkl?

Re:Pardiso got the wrong Schur complement

Mahan — Wed, 04 Sep 2024 05:02:57 GMT

Hello Brian

Could you please calculate and provide the condition number for this matrix?

Re: Re:Pardiso got the wrong Schur complement

BrianZHANG — Wed, 04 Sep 2024 08:11:39 GMT

Hi Mahan

I calculated the condition number of this matrix using Matlab's condest, and the result is 1.52*10^14. I know that the ill-conditioned system can lead to numerical errors, but why are the non-zero patterns of the two so different, especially in the upper right corner of the matrix? The result of Pardiso looks like some columns are shifted to the left.

When I do Jacobi preprocessing, which is multiplying by a diagonal matrix, the condition number is reduced to 90,000. I have set iparm(11)=1, iparm(13)=1, and iparm(36)=2 to compute dense schur's complement, and the result is still the same. Scaling and weighted matching should also help improve matrix properties, shouldn't they?

Re:Pardiso got the wrong Schur complement

Mahan — Fri, 13 Sep 2024 02:13:47 GMT

Hi Brian,

Could you please share a reproducer for this? The example Fortran source you have shared is not compilable.

Re: Re:Pardiso got the wrong Schur complement

BrianZHANG — Fri, 13 Sep 2024 07:05:17 GMT

Hi Mahan,

Thank you for your time to check it for me. I uploaded the VS project and Matlab code. The original matrix was too large, so I uploaded it to Google Drive. The matrix after jacobi preprocessing is also included.

Re:Pardiso got the wrong Schur complement

Mahan — Fri, 27 Sep 2024 06:50:29 GMT

Hi Brian,

The error is reproducible and the issue has been communicated to the development team. Please allow us a few days to come up with an update on this.

Re:Pardiso got the wrong Schur complement

Mahan — Fri, 04 Oct 2024 02:14:06 GMT

Hi Brian,

This issue might take some time to resolve, I am not having any specific ETA for this at this moment.

Re: Re:Pardiso got the wrong Schur complement

BrianZHANG — Wed, 16 Oct 2024 02:13:38 GMT

Hi Mahan,
I'll await further updates and appreciate any information when available. Thanks again for your support.

Re: Pardiso got the wrong Schur complement

Chao_Y_Intel — Fri, 06 Feb 2026 04:47:05 GMT

Brain, can you please check with new MKL releases: https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html This problem expected to be fixed starting from MKL 2025.2. Thanks for your report of this problem.

thanks,
Chao