MKL/Pardiso sensitive to zero-valued entries (ZVE)

[mecej4]

When MKL/Pardiso is used to solve a system of N X N sparse equations, it may happen that some of the sparse matrix entries have zero values. The presence of such zero values may be inadvertent or the result of a finite-element assembly process, or they may be present in a standard test matrix obtained from a repository.  When a large number of such zero-valued entries (ZVE) are present, it is understandable that the bandwidth and solution may be affected, but it is not expected that the solution should be affected by the presence of a handful of ZVEs except for small changes attributable to finite-precision arithmetic.

Here is a problem reproducer, based on the matrix IMPCOL_C from the NIST Matrix Market ( http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/chemimp/impcol_c.html ). This unsymmetric matrix is 137 X 137, with 411 entries, of which 11 are ZVEs. For the purposes of this report, I removed 8 of these 11, plus another non-zero, leaving only three ZVEs in the matrix (data and test program file source attached). The test program reads the data in CSR format, and checks the solution against the known vector x₀ that was used to obtain the r.h.s. vector b = Ax₀. Likewise, the program checks the residuals Ax - b and prints those above a certain threshold.

With the three ZVEs kept in the matrix, the program output is:

s:\DARIO\UNSYMM>unsymmkl impcolcA.csr
Intel(R) Math Kernel Library Version 11.2.4 Product Build 20150805 for 32-bit applications

Matrix is 137 X 137, with 402 entries (of which 3 are zero)
Location of zero-valued entries:
( 51, 58)  ( 81, 88)  ( 84, 91)
Pardiso iparm[] non-zero entries:
  p[ 0]= 1  p[ 1]= 2  p[ 2]= 1  p[ 9]=13  p[10]= 1  p[12]= 1  p[17]=-1  p[18]=-1


Reordering completed ...
Number of nonzeros in factors  = 1017
Factorization completed ...
Solve completed ...

Solution components with error > 1.00e-06, if any:
  i     x(i)       x(i)/x_0(i) - 1
  24   1.16912e+00  -7.265e-06
 115   1.83824e+00  -7.265e-06
 122   1.88971e+00  -4.039e-06

Residuals, if any > 1d-4
    .... None
Maximum Residual =  7.872e-06
Norm of Residual =  1.124e-05

The same program/problem combination, rerun with the -z option to cause removal of the three ZVEs before calling Pardiso, gives:

s:\DARIO\UNSYMM>unsymmkl impcolcA.csr -z
Intel(R) Math Kernel Library Version 11.2.4 Product Build 20150805 for 32-bit applications
                                                            
Matrix is 137 X 137, with 402 entries (of which 3 are zero)
Location of zero-valued entries:
( 51, 58)  ( 81, 88)  ( 84, 91)
RMZEROS: nnz reduced from 402 to 399
Pardiso iparm[] non-zero entries:
  p[ 0]= 1  p[ 1]= 2  p[ 2]= 1  p[ 9]=13  p[10]= 1  p[12]= 1  p[17]=-1  p[18]=-1


Reordering completed ...
Number of nonzeros in factors  = 1185
Factorization completed ...
Solve completed ...

Solution components with error > 1.00e-006, if any:
    .... None

Residuals, if any > 1d-4
    .... None
Maximum Residual = 3.553e-015
Norm of Residual = 1.439e-014

Note the drastic change in the residual norm. In fact, user Dario Mangoni gave in https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/592612 some example matrices for which the residuals changed from 10⁹ to 10^-15 as a resulting of removing ZVEs before solving.

Similar results are obtained if the matrix IMPCOL_C is used as-is (i.e., with the 11 ZVEs intact).

I have found that the third party libraries Pardiso-4.12 and Pardiso-5 from Basel/Lugano have this problem, as well. At the least, Intel/MKL should caution users to remove ZVEs from their matrix data before calling Pardiso. It would be nice to provide a utility function to do this in MKL, or you could make ZVE removal part of the default behavior. In the long run, however, the unexpected sensitivity of Pardiso to the presence of ZVEs should be investigated and fixed.

The matrix in question is somewhat ill-conditioned, as were Dario Mangoni's. However, the condition number of IMPCOL_C is about 10⁵, so the condition number does not seem solely responsible for the inaccurate solutions and large residuals.

Thank you.

[Alexander_K_Intel2]

Hi Mecej4,

Thanks a lot for this interesting analysis. However I don't think that in general removing zeros elements from matrix helps to get stable result. It can either increase or decrease overall residual. Removing zero elements change portrait of matrix and it's permutation array. We cannot operate with theorem (||r||<condA*||f||) because it take place for algorithms with global pivoting which is impossible for multifrontal methods. 

Thanks,

Alex