Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Every other call to dss_solve_complex gives incorrect answer

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

Paul_N_2

Beginner

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

04-16-2016
08:55 PM

59 Views

Every other call to dss_solve_complex gives incorrect answer

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

5 Replies

mecej4

Black Belt

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

04-17-2016
07:18 PM

59 Views

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

Paul_N_2

Beginner

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

04-17-2016
07:35 PM

59 Views

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

Paul_N_2

Beginner

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

04-17-2016
08:25 PM

59 Views

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

mecej4

Black Belt

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

04-17-2016
08:29 PM

59 Views

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.

mecej4

Black Belt

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

04-17-2016
08:32 PM

59 Views

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

For more complete information about compiler optimizations, see our Optimization Notice.