Pardiso got the wrong Schur complement

[BrianZHANG]

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?

[Mahan]

Hello Brian

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

[BrianZHANG]

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?

[Mahan]

Hi Brian,

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

[BrianZHANG]

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.

[Mahan]

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.

[Mahan]

Hi Brian,

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

[BrianZHANG]

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

[Chao_Y_Intel]

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