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?
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
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.
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.
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, to the program. The file should contain compressed sparse row matrix data in the following format, with all indices 1-based.
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.
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
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?
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 gives the number of positive eigenvalues and and iparm 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.
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)
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
1-st repetition. Inertia : 73 6 0
2-nd repetition. Inertia : 73 6 0
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
1-st repetition. Inertia : 5 3 0
2-nd repetition. Inertia : 3 5 0
3-rd repetition. Inertia : 3 5 0
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