- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I am trying to call the DSS solver in MKL from the Intel Python distribution 2017 beta via the ctypes library. I have everything working properly except for repeated calls to dss_solve_complex. The first, third, fifth, ... calls yields the correct solution, yet the second, forth,... and so on give the wrong answer. Moreoever, the output is not the same for repeated runs. In addition, the solver returns success as the output even though the answer is incorrect. This same issue does not occur when calling dss_solve_real. A minimal example that can be used in both the real and complex solvers is below:

A

===

[[ 1 0 2]

[ 0 0 3]

[-4 5 6]]

b

===

[0 2 0]

The correct solution should be: [-1.33333333 -1.86666667 0.66666667]

Trying to repeatedly solve five times yields, for example:

[-1.33333333+0.j -1.86666667+0.j 0.66666667+0.j]

[ 3.13043594e-04 +3.08799744e-03j 1.50203705e-05 +1.60000000e+01j

1.60000000e+01 +1.00000000e+00j]

[-1.33333333+0.j -1.86666667+0.j 0.66666667+0.j]

[ 3.13043594e-04 +3.08799744e-03j 1.19209290e-05 +1.60000000e+01j

1.60000000e+01 +1.00000000e+00j]

[-1.33333333+0.j -1.86666667+0.j 0.66666667+0.j]

Again, I have tried to run this in the Intel Python Beta 2017, and also the latest Anaconda Python distro that uses MKL 11.3.1. The Python code that is used to call the DSS routines can be found here:

https://github.com/nonhermitian/qutip/blob/mkl_enhancements/qutip/mkl/spsolve.py

Best,

Paul Nation

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Such incorrect results can originate in (i) bugs in user code (ii) bugs in non-native language interface code (iii) bugs in the MKL library routines. I am not a Python user, nor am I up to looking at nearly 400 lines of Python code to find any errors.

However, a short Fortran program modeled after the example code provided with MKL enables me to discount (iii). The program works on your test matrix, repeating seven cycles of three different sets of right hand side vectors. The r.h.s. vectors for irep = 1, 4, 7, ..., 19 are the same as the one that you specified.

! ifort /Qmkl /nostand /warn:none dss_cmplx.f90 PROGRAM dss_test USE mkl_dss IMPLICIT NONE INTEGER, PARAMETER :: nrows = 3, ncols = 3, nnonzeros = 7, nrhs = 1 INTEGER :: rowindex(nrows+1) = [ 1, 3, 5, 8 ] INTEGER :: columns(nnonzeros) = [ 1, 3, 2, 3, 1, 2, 3 ] DOUBLE COMPLEX :: values(nnonzeros) = [ (1.,0.), (2.,0.), (0.,0.), & (3.,0.), (-4.,0.), (5.,0.), (6.,0.) ], rhs0(nrows) = [ (0.,0.), (2., & 0.), (0.,0.) ], rhs(nrows) DOUBLE COMPLEX solution(nrows) INTEGER *8 handle INTEGER i, error, buf, idum(1), irep error = dss_create(handle, mkl_dss_defaults) IF (error/=mkl_dss_success) GO TO 100 error = dss_define_structure(handle, mkl_dss_non_symmetric, rowindex, & nrows, ncols, columns, nnonzeros) IF (error/=mkl_dss_success) GO TO 100 error = dss_reorder(handle, mkl_dss_defaults, idum) IF (error/=mkl_dss_success) GO TO 100 error = dss_factor_complex(handle, mkl_dss_indefinite, values) IF (error/=mkl_dss_success) GO TO 100 DO irep = 1, 21 rhs(1:3) = rhs0(1:3) + mod(irep-1, 3)* [ (0.3,0.0), (0.2,0.0), & (-0.2, 0.) ] error = dss_solve_complex(handle, mkl_dss_defaults, rhs, nrhs, & solution) IF (error/=mkl_dss_success) GO TO 100 IF (mod(irep,3)==1) WRITE (*, *) WRITE (*, '(i2,3(2x,(F8.4,F8.4)))') irep, (solution(i), i=1, ncols) END DO error = dss_delete(handle, mkl_dss_defaults) IF (error/=mkl_dss_success) GO TO 100 STOP 100 WRITE (*, *) 'DSS error code ', error STOP 1 END PROGRAM dss_test

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Many thanks for taking a look at this. I guess there is some issue with my call somewhere. One question though, why can you use MKL_DSS_NON_SYMMETRIC rather than MKL_DSS_NON_SYMMETRIC_COMPLEX in the call to define_structure? Isn't that for real problems only?

Best,

Paul

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Ok, I found the error. It turns out that asking for the statistics after each call to solve was some how messing up the next call. Now everything works great! Thanks again for taking a look at this. It helped a lot.

Best,

Paul

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Good question!

I was not aware that definitions with _COMPLEX at the end had been added! DSS has been in existence for over 15 years, and I had adapted code that came with Compaq/Digital Fortran (circa 2000). In the CVF manual there is no mention of MKL_DSS_SYMMETRIC_STRUCTURE_COMPLEX.

Since the DSS in MKL is just a wrapper around Pardiso, I find it unsurprising that my error (writing MKL_DSS_SYMMETRIC_STRUCTURE instead) did not cause problems. Pardiso uses the CSR representation of sparse matrices, and the matrix is described in terms of indices rather than C-pointers. Therefore, the size of each matrix element should have no effect on the CSR data.

We may hope that one of the Intel-MKL personnel will clarify this point.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Paul N. wrote:

It turns out that asking for the statistics after each call to solve was some how messing up the next call. Now everything works great!

That's actually bad news, because of the implication that the bug can go into hiding as a result of removing the statistics inquiry calls. It may be worth finding out if the same thing happens with the Fortran code.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page