Can I multiply the real parts of two complex matrices using dgemm?

Question

I have two complex matrices A and B, with matching shapes.

Is there a way to cleverly setup the dgemm arguments so as to get the result of the matrix multiplications of the real parts of these matrices only?

Intuitively, the real and imaginary parts of a complex matrix are stored in an interlaced way, so I guess they correspond to strided real arrays.

Note: I am using cython and accessing blas via scipy.linalg.cython_blas, but any pure C solution is welcome as well.

Answer 1

After some research, I think this is impossible using standard BLAS, but I was able to do this using BLIS. Here is an example in Cython:

cimport blis.cy
from libc.stdlib cimport abort, malloc, free

# Computes the real part of a @ b
cpdef void blis_matmul_real(
    double complex[:,::1] a ,
    double complex[:,::1] b ,
    double[:,::1] c         ,
) noexcept nogil:

    cdef int m = a.shape[0]
    cdef int k = a.shape[1]
    cdef int n = b.shape[1]

    cdef double* a_real = <double*> &a[0,0]
    cdef double* b_real = <double*> &b[0,0]
    cdef double* res = &c[0,0]

    cdef int lda = 2 * k
    cdef int ldb = 2 * n

    blis.cy.gemm(
        blis.cy.NO_TRANSPOSE, blis.cy.NO_TRANSPOSE,
        m, n, k,
        1.0, a_real, lda, 2,
        b_real, ldb, 2 ,
        0.0, res, n, 1
    )

    # Pointer addition
    a_real += 1
    b_real += 1

    blis.cy.gemm(
        blis.cy.NO_TRANSPOSE, blis.cy.NO_TRANSPOSE,
        m, n, k,
        -1.0, a_real, lda, 2,
        b_real, ldb, 2 ,
        1.0, res, n, 1
    )