I'm using Pardiso to compute the Schur complement of a symmetric positive-semidefinite matrix with a kernel of dimension 6. I don't know if this is officially supported, but for similar rank-deficient matrices, I sometimes get the correct results. To see whether the Schur complement is right or not, I'm checking in MATLAB the eigenvalues of the result given by Pardiso after the numerical factorization. The first 6 eigenvalues should be 0. With MUMPS, I get the following (correct) output:
>> eig(triu(MUMPS) + tril(MUMPS') - diag(diag(MUMPS)))
With brute force MATLAB, I get the following (correct) output:
>> sort(real(eig(MATLAB(10:63,10:63) - MATLAB(10:63,1:9) * MATLAB(1:9,1:9)^-1 * MATLAB(1:9,10:63))))
With Pardiso, I get the following (wrong) output:
>> eig(triu(Pardiso) + tril(Pardiso') - diag(diag(Pardiso)))
I'm attaching a reproducer that factorizes the matrix and output the Schur complement, as well as the original matrix stored in dense format to be used inside MATLAB. Can you reproduce this problem ? Do I need to adjust some parameters, or is there a bug inside the solver ?
Thank you for your help.
Your test returns -1, which means the input inconsistency. I enabled matrix checker and it reported more details, see below
*** Error in PARDISO ( sequence_ido,parameters) error_num= 10
*** Input check: iparam -858993460 (out of bounds)
*** Input parameters: inconsistent error= 10 max_fac_store_in: 1
matrix_number_in: 1 matrix_type_in: -2
ido_in : 12 neqns_in : 63
ia[neqns_in]-1 : 707 nb_in : 1
I have a problem with Schur complement. In some cases when I calculate Schur complement I get wrong matrix. Values in matrix are correct but on the wrong place in the matrix. Error is appearing randomly.
Any ideas are welcome, thanks.
Nina, could you try the latest update 2 ( MKL 2017)? We have fixed the issue when MKL Pardiso Schur complement returns incorrect results. see MKL v.2017 Bug Fix list - https://software.intel.com/en-us/articles/intel-math-kernel-library-intel-mkl-2017-bug-fixes-list
Nina Z. wrote:
With latest update it works now for symmetric matrix, but for unsymmetrical there is the same bug as before for symmetric.
Can you please provide us with a reproducer for nonsymmetric case?