Complex inner products always return zero

[jmbr]

The outcome of

(magicl.blas-cffi::%ddot 1 (make-array 1 :element-type 'double-float :initial-element 1.0d0) 1 (make-array 1 :element-type 'double-float :initial-element -1.0d0) 1)

is -1.0d0, as expected. However, the outcome of

(magicl.blas-cffi::%zdotu 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(1.0d0 0.0d0)) 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(-1.0d0 0.0d0)) 1)

is #C(0.0d0 0.0d0). The resulting zero value is consistent no matter which arrays are passed as input.

[jmbr]

The problem stems from the fact that there exist two types of application binary interfaces (ABIs) for BLAS:

1. On the one hand, the reference implementation of BLAS as well as OpenBLAS both follow the standard C ABI [1] and return complex values via the SSE registers XMM0 (real part) and XMM1 (imaginary part). When using those implementations, MAGICL returns the correct dot product.
2. On the other hand, Apple's Accelerate framework and the Intel Math Kernel library deviate from the standard [2] and modify the interface so that functions that should return a complex number instead have a void return type and an extra argument that points to the memory where the complex result will be written.

To summarize: MAGICL's interface to zdotu, cdotu, etc. is fine if the underlying BLAS is OpenBLAS or the reference BLAS but things fail as shown above when using Apple Accelerate or Intel MKL. In the latter cases, we need to define a different interface to make things work:

(CFFI:DEFCFUN ("zdotu_" %%ZDOTU :LIBRARY MAGICL.FOREIGN-LIBRARIES:LIBBLAS)
    :VOID
  (C :POINTER)
  (N :POINTER)
  (ZX :POINTER)
  (INCX :POINTER)
  (ZY :POINTER)
  (INCY :POINTER))

(COMMON-LISP:DEFUN %ZDOTU (N ZX INCX ZY INCY)
  (COMMON-LISP:DECLARE (COMMON-LISP:INLINE %%ZDOTU)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) N)
                       (COMMON-LISP:TYPE
                        (COMMON-LISP:SIMPLE-ARRAY
                         (COMMON-LISP:COMPLEX COMMON-LISP:DOUBLE-FLOAT)
                         (COMMON-LISP:*))
                        ZX)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) INCX)
                       (COMMON-LISP:TYPE
                        (COMMON-LISP:SIMPLE-ARRAY
                         (COMMON-LISP:COMPLEX COMMON-LISP:DOUBLE-FLOAT)
                         (COMMON-LISP:*))
                        ZY)
                       (COMMON-LISP:TYPE (COMMON-LISP:SIGNED-BYTE 32) INCY))
  (CFFI:WITH-FOREIGN-OBJECTS ((PTR-C :DOUBLE 2) (N-REF885 ':INT32) (INCX-REF887 ':INT32)
                              (INCY-REF889 ':INT32))
    (COMMON-LISP:SETF (CFFI:MEM-REF N-REF885 :INT32) N)
    (COMMON-LISP:SETF (CFFI:MEM-REF INCX-REF887 :INT32) INCX)
    (COMMON-LISP:SETF (CFFI:MEM-REF INCY-REF889 :INT32) INCY)
    (MAGICL.CFFI-TYPES:WITH-ARRAY-POINTERS ((ZX-REF886 ZX) (ZY-REF888 ZY))
      (%%ZDOTU PTR-C N-REF885 ZX-REF886 INCX-REF887 ZY-REF888 INCY-REF889)
      (COMMON-LISP:COMPLEX (CFFI:MEM-AREF PTR-C :DOUBLE 0)
                           (CFFI:MEM-AREF PTR-C :DOUBLE 1)))))

[1] https://github.com/hjl-tools/x86-psABI/wiki/X86-psABI [2] https://software.intel.com/en-us/mkl-macos-developer-guide-calling-blas-functions-that-return-the-complex-values-in-c-c-code

[appleby]

Nice sleuthing.

[jmbr]

Thank you.

[colescott]

Alright I was able to reproduce this with MKL and got a fun result:

* (magicl.blas-cffi::%ddot 1 (make-array 1 :element-type 'double-float :initial-element 1.0d0) 1 (make-array 1 :element-type 'double-float :initial-element -1.0d0) 1)
-1.0d0
* (magicl.blas-cffi::%zdotu 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(1.0d0 0.0d0)) 1 (make-array 1 :element-type '(complex double-float) :initial-element #c(-1.0d0 0.0d0)) 1)
CORRUPTION WARNING in SBCL pid 8901(tid 0x7fc6591390c0):
Memory fault at (nil) (pc=0x7fc643cc728a, fp=0x7fc64bcff3f0, sp=0x7fc64bcff2c0) tid 0x7fc6591390c0
The integrity of this image is possibly compromised.
Continuing with fingers crossed.

debugger invoked on a SB-SYS:MEMORY-FAULT-ERROR in thread
#<THREAD "main thread" RUNNING {10005184C3}>:
  Unhandled memory fault at #x0.

I'll go ahead and try regenerating the bindings and give an update here.