MKL DSS results change in 17.0.1

[AndrewC]

I have upgraded from Intel Parallel Studio 15.0 to 17.0 Update 1 and am noticing that MKL DSS (aka PARDISO) is giving somewhat different results when solving large real, symmetric sparse matrices. Is this a known issue? I am seeing the differences on both Windows and Linux. For example, has the reordering strategy changed or something like that?

[AndrewC]

The matrices are 79x79. That's what I am passing to MKL.  And yes, the last row is all 0. Since the matrix is small it should be easy to verify the ROW COMPRESSED format being passed to DSS by looking at the p and i arrays.

[EDIT] Actually, sorry about that. the above statement is incorrect. The last row is not all 0. If you look at the matrix file , you will see an entry 78,78 below which is the diagonal entry. There is no requirement that the i,j,values be in increasing order of i,j

 78 78 0.042913080130498

[mecej4]

I see. I don't know what the csparse package does in detail, but as long as you know that it reads the file and puts out the CSR data correctly, that should be OK. 

[AndrewC]

I've been using CSparse for a while and it's solid.

I also loaded these matrices into Matlab and confirmed the # of eigenvalues. In particular there are 6 -ve eigenvalues for matrix #2.

[mecej4]

It seems to me that some memory corruption is happening in dss_statistics(). In your routine full_factor_matrix(), if statIn is set to "Inertia,FactorTime", the second next call to dss_statistics(), which is in main(), returns the correct values for the inertia. If we add ",flops" to statIn in full_factor_matrix(), the second next call to dss_statistics() returns incorrect values. In both cases, the first call to dss_statistics() returns the correct results.

If, as I think, dss_statistics() merely reports results that are already available, selecting the requested items from the whole set, it is reasonable to suspect that memory corruption occurred in the first call to dss_statistics().

Since the bug is not always seen (on a different PC with VS 2015, I did not see it), it may be a bit tough to track the source of the error.

As a workaround, you could specify only "Inertia", if that will suffice and if this change fixes not only the test code but also your real application.

[mecej4]

I have created a simplified reproducer for the problem which is self-contained (does not need the CSPARSE package) and exhibits the bug even with small test matrices. The bug is present in MKL 11.3.4 and 2017.0.1, 32- and 64-bit, on Windows 10 Pro 64-bit, with VS2013. The bug is not present in MKL 11.1.4 and 11.2.4 on the same machine.

The issue is related to factorizing a symmetric indefinite matrix and obtaining the inertia of the matrix by making calls to the DSS sparse solver using the recommended procedure. The last of these calls is a call to dss_statistics(). If, in this call, the request in statin is just "Inertia", no bug occurs at all. If the request contains something more, for example, "Inertia,Flops", the first call to dss_statistics() returns the correct inertia, but subsequent pairs of calls to dss_factor() and dss_statistics(), with different matrix entry values, or even the same matrix values, return incorrect results.

The reproducer reads a test matrix from a file whose name is given as the first argument, argv[1], to the program. The file should contain compressed sparse row matrix data in the following format, with all indices 1-based.

1. Line-1: n, nnz (matrix size is n X n, matrix has nnz non-zero entries)
2. Lines 2 to n+2: row indices in CSR representation
3. Lines n+2 to n+1+nnz: (column, value) pairs for non-zero entries

Three test matrices are provided: s.csr contains data for a small 8 X 8 matrix with 18 non-zero entries. The files mat1.csr and mat2.csr contain the same matrices as in vascii_'s Matijv_1.txt and Matijv_2.txt, but in CSR format and with 1-based indices. I suspect that any other symmetric indefinite matrix will also serve just as well.

After compiling and linking the source dsscsr.c, run the examples using the batch file run.bat (or otherwise, as you prefer). Here is a typical output:

s:\lang\MKL\DSSTST>dsscsr mat2.csr
Intel(R) Math Kernel Library Version 2017.0.1 Product Build 20161005 for 32-bit applications
Matrix has 79 rows, 79 columns, 425 non-zero entries
Calling DSS_CREATE
Calling DSS_DEFINE
Calling DSS_REORDER
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  1-st repetition. Inertia :    73     6     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  2-nd repetition. Inertia :    65    14     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  3-rd repetition. Inertia :    65    14     0

Only the first result (73,6,0) is correct. The other results are the same incorrect results that vascii_ observed.

Change "Inertia,Flops" to "Inertia" in the source file, rebuild and rerun to see that this change causes the bug to go away.

The bug is a bit elusive. On a Windows 10 Home system, with VS2015, the bug occurs only with the small matrix in s.csr and not with the other two matrix files.

[mecej4]

Here is a Fortran version of the test program of #26. On both the systems that I used (one has VS2013, the other VS2015; both have various versions of MKL), I did not see the bug that the C version showed.

Compile with MKL and, one at a time, redirect the CSR data files provided in #26 to the standard input of the program.

! Use DSS to solve N X N linear equations with solution 1..2.
! Matrix is assumed to be symmetric indefinite and only the 
! upper triangle is input from a CSR data file containing
!
! Line 1: n,nnz
! Lines 2 to (n+2): row indices, 1-based
! Lines n+3 to n+2+nnz: column, value
!       (column indices 1-based, lines in required order)
!
    program dss_test
      use mkl_dss
      implicit none
!
      integer :: nrow, ncol, nnz
      integer, allocatable :: rowindex(:), columns(:)
      double precision, allocatable :: values(:)
      type (mkl_dss_handle) :: handle
      integer error
      integer i, j, k, irep, idum(1)
      logical, allocatable :: empty(:)
      double precision :: statout(4)
      character*40 :: statin
      integer buf(10)
      equivalence (buf, statin)
!	  
      read (*, *) nrow, nnz
      allocate (rowindex(nrow+1), empty(nrow))
      read (*, *) rowindex(1:nrow+1)
      ncol = nrow
      allocate (columns(nnz), values(nnz))
      read (*, *)(columns(i), values(i), i=1, nnz)
!	  rowindex(1:nrow+1) = rowindex(1:nrow+1)+1 ! 0-base to 1-base
!	  columns(1:nnz)=columns(1:nnz)+1           ! 0-base to 1-base
      empty = .true.
      do k = 1, nnz
        j = columns(k)
        if (values(k)/=0d0) empty(j) = .false.
      end do
      do i = 1, nrow
        if (empty(i)) write (*, '(1x,A,1x,i6,1x,A)') 'Col', i, 'is empty'
      end do
      write (*, 100) 'Matrix has ', nrow, ' rows ', nrow, nnz
100   format (1x, a, i5, a, i5, ' columns and ', i5, ' non-zero entries')
!
      print *, ' calling DSS_CREATE'
      error = dss_create(handle, mkl_dss_defaults)
      if (error/=mkl_dss_success) go to 110
!
      print *, ' calling DSS_DEFINE'
      error = dss_define_structure(handle, mkl_dss_symmetric, rowindex, nrow, &
        ncol, columns, nnz)
      if (error/=mkl_dss_success) go to 110
      print *, ' calling DSS_REORDER'
      error = dss_reorder(handle, mkl_dss_defaults, idum)
      if (error/=mkl_dss_success) go to 110
      print *, ' calling DSS_FACTOR'
      error = dss_factor_real(handle, mkl_dss_indefinite, values)
      if (error/=mkl_dss_success) go to 110
      do irep = 1, 3
        print *, ' calling DSS_STATISTICS'
        statin = 'INERTIA,flops'C
        error = dss_statistics(handle, mkl_dss_defaults, buf, statout)
        if (error/=mkl_dss_success) go to 110
        write (*, '(10x,A,2x,3F6.0)') 'Inertia : ', statout(1:3)
      end do
      error = dss_delete(handle, mkl_dss_defaults)
      if (error/=mkl_dss_success) go to 110
      stop
110   write (*, *) 'DSS error code ', error
      stop 1
    end program dss_test

[AndrewC]

Great stuff , thanks for going to all the work to create the various reproducers!

I'll check out the suggested work-around of only asking for 'Inertia'. This brings back some old memories of problems with dss_statistics in older versions of MKL. Perhaps MKL 10.0? Asking for more than 'Interia' would cause a crash or return rubbish values. This may be an issue that has been bouncing around waiting to bite.

Sounds like some kind of subtle memory overwrite error. Perhaps passing C strings to FORTRAN?
 

[AndrewC]

Update: changing the string to "Inertia" ONLY seems to correct the problems I am seeing in my QA cases.

[mecej4]

Another approach for you to consider using is to call Pardiso directly, rather than through the DSS wrappers. After setting all the parameters for Pardiso, you only need to call Pardiso twice, with phase=11 once, and then repeatedly call with 22 for each set of matrix values with the same structure. After the phase=22 call, iparm[21]  gives the number of positive eigenvalues and and iparm[22] contains number of negative eigenvalues.

The DSS existed in the DEC/Compaq DXML/CXML library (and still exists in the Sun/Oracle Performance  Library), but it is now more than a decade since marketing of CXML was stopped.

I realize that if your code contains DSS calls in several places, you may not want to do the conversion.

[Gennady_F_Intel]

Colleagues, that's really elusive case. 

I tried 32 bit MKL 2017 with mat2.csr and s.csr inputs. 

with mat2.csr, I see on my side ( win 8.1, HSW, 2 cores)

C:\5_Customers\ESI\ww01_2017_#1>32.exe mat2.csr
Intel(R) Math Kernel Library Version 2017.0.1 Product Build 20161005 for 32-bit applications
Matrix has 79 rows, 79 columns, 425 non-zero entries
Calling DSS_CREATE
Calling DSS_DEFINE
Calling DSS_REORDER
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  1-st repetition. Inertia :    73     6     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  2-nd repetition. Inertia :    73     6     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  3-rd repetition. Inertia :    73     6     0

but with s.csr I see the follow result which I believe is wrong:

Intel(R) Math Kernel Library Version 2017.0.1 Product Build 20161005 for 32-bit applications
Matrix has 8 rows, 8 columns, 18 non-zero entries
Calling DSS_CREATE
Calling DSS_DEFINE
Calling DSS_REORDER
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  1-st repetition. Inertia :     5     3     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  2-nd repetition. Inertia :     3     5     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  3-rd repetition. Inertia :     3     5     0

wbr, Gennady

[mecej4]

Gennady, I am glad that you can reproduce at least one instance of the bug! The data in s.csr is taken from Olaf Schenk's example for Pardiso, http://www.pardiso-project.org/manual/pardiso_sym.f . The Pardiso site has a C version, as well, with 0-based indices.

I was slightly desperate after trying to reproduce the problem with several medium sized symmetric indefinite matrices from the Matrix Market, and failing. The small matrix offered a case where typing in the data could be done without spending a lot of time.

The following is an adaptation of one of the matrix generators from the Matrix Market. It takes two arguments: the number of rows (= number of columns) of the matrix, and the sparseness (real number: 0 for dense matrix, 1 for diagonal matrix). There is no guarantee that the generated matrices will have at least some negative eigenvalues, but it seemed to work fine for several inputs that I tried.

! Purpose: generate a symmetric matrix file in Compressed Sparse Row 
!          for use with DSSCSR tests.
!
    Program dmatgen

      Implicit None

      Integer nmax
      Parameter (nmax=500)

      Character :: dist = 'N', rsign = 'F', sym = 'S', &
                   grade = 'N', pivtng = 'N', pack = 'R'
      Integer :: info, kl = 9, ku = 9, lda, m, n
      Integer :: mode = 6, model = 6, moder = 6
      Integer :: nnz, nnzreq, i, j, idx
      Double Precision :: cond = 1, condr = 1, condl = 1.0, anorm, sparse
      Double Precision dmax, v

      Integer iseedout(4), iseedin(4), ipivot(nmax)
      Integer iwork(nmax), ia(nmax+1), col(nmax*(nmax+1)/2)
      Double Precision a(nmax*nmax)
      Double Precision d(nmax), dl(nmax), dr(nmax)

      iseedin = [2264, 1462, 3476, 3231]
      Write (*, '(A,$)') 'Enter m, sparseness : '
      Read (*, *) m, sparse
      Open (Unit=10, File='sym.csr', Status='replace')
      n = m
      dmax = 1
      anorm = -1
      lda = m
!
!     'iseedin' is used to save a copy of the seed values for printout
!     along with the matrix data, since dlatmr overwrites the
!     seed argument with values to be used on the next call.
!
      iseedout(1:4) = iseedin(1:4)

      Call dlatmr(m, n, dist, iseedout, sym, d, mode, cond, dmax, rsign, &
        grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, &
        sparse, anorm, pack, a, lda, iwork, info)

      If (info/=0) Then

        Print *, 'Error in subroutine dlatmr: ', info
        Stop

      End If

      nnzreq = (m*n-m)/2 + m
!
!       Convert from lower triangle to CSR upper triangle output format:
!
      nnz = 0
      idx = 0
      Do j = 1, n
        idx = idx + 1
        v = a(idx)
        If (v==0) v = 1.5 + 0.3*abs(2*j-n)/dble(n)
        nnz = nnz + 1
        ia(j) = nnz
        a(nnz) = v
        col(nnz) = j
        Do i = j + 1, n
          idx = idx + 1
          If (a(idx)/=0.0) Then
            nnz = nnz + 1
            a(nnz) = a(idx)
            col(nnz) = i
          End If
        End Do
      End Do
      ia(n+1) = nnz + 1
      Write (10, '(1x,2I8)') n, nnz
      Write (10, '(1x,I8)')(ia(i), i=1, n+1)
      Write (10, '(1x,I8,1x,1pE13.5)')(col(i), a(i), i=1, nnz)
      Close (10)
      Stop
    End Program dmatgen

As I reported earlier, the three matrices in s.csr, mat1.csr and mat2.csr did not result in the bug when I compiled under VS2015. Using this generator, I am now able to turn up a few cases. Run the generator with m = 8, sparseness = 0.5, or with m = 9, 0.5. For the last pair, DSS reports:

S:\LANG\MKL>dsymgen & dsscsr sym.csr
Enter m, sparseness : 9,0.5
Intel(R) Math Kernel Library Version 11.3.4 Product Build 20160823 for 32-bit applications
Matrix has 9 rows, 9 columns, 25 non-zero entries
Calling DSS_CREATE
Calling DSS_DEFINE
Calling DSS_REORDER
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  1-st repetition. Inertia :     7     2     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  2-nd repetition. Inertia :     6     3     0
Calling DSS_FACTOR_REAL
Calling DSS_STATISTICS
  3-rd repetition. Inertia :     6     3     0