Every other call to dss_solve_complex gives incorrect answer

[Paul_N_2]

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

[mecej4]

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]

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]

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]

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]

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.