能否通过配置dgemm参数实现两个同形复矩阵实部的矩阵乘法？

仅通过配置dgemm的参数无法直接实现该需求。BLAS的dgemm仅支持处理连续存储的实矩阵块，不提供元素级别的步长设置——而复矩阵的实部元素在内存中是与虚部交错存储的（每个复数占两个连续double，实部在前、虚部在后），无法被dgemm直接识别为连续的矩阵结构。

不过可以通过以下两种可行方案实现目标：

方案1：提取实部到连续矩阵后调用dgemm #

这是最直接、易维护的方案，内存拷贝的开销通常远小于矩阵乘法的计算开销。

C语言示例（使用CBLAS） #

#include <complex.h>
#include <cblas.h>
#include <stdlib.h>

void multiply_complex_realparts(int m, int k, int n, 
                                double complex *A, double complex *B, double *C) {
    // 分配实部矩阵内存（列优先存储）
    double *A_real = malloc(m * k * sizeof(double));
    double *B_real = malloc(k * n * sizeof(double));
    
    // 提取A的实部
    for (int j = 0; j < k; j++) {
        for (int i = 0; i < m; i++) {
            A_real[i + j * m] = creal(A[i + j * m]);
        }
    }
    
    // 提取B的实部
    for (int j = 0; j < n; j++) {
        for (int t = 0; t < k; t++) {
            B_real[t + j * k] = creal(B[t + j * k]);
        }
    }
    
    // 调用dgemm计算实部矩阵乘法
    const char transa = 'N', transb = 'N';
    const double alpha = 1.0, beta = 0.0;
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                m, n, k, alpha, A_real, m, B_real, k, beta, C, m);
    
    // 释放临时内存
    free(A_real);
    free(B_real);
}

Cython示例（使用scipy.linalg.cython_blas） #

from scipy.linalg.cython_blas cimport dgemm
import numpy as np
cimport numpy as np

def multiply_complex_realparts(np.ndarray[np.complex128_t, ndim=2] A, 
                               np.ndarray[np.complex128_t, ndim=2] B):
    cdef int m = A.shape[0], k = A.shape[1], n = B.shape[1]
    cdef np.ndarray[np.float64_t, ndim=2] A_real = A.real.astype(np.float64)
    cdef np.ndarray[np.float64_t, ndim=2] B_real = B.real.astype(np.float64)
    cdef np.ndarray[np.float64_t, ndim=2] C = np.zeros((m, n), dtype=np.float64)
    
    # 调用dgemm（遵循CBLAS列优先格式）
    dgemm(b'N', b'N', &m, &n, &k, 1.0, 
          <double*>A_real.data, &m, 
          <double*>B_real.data, &k, 
          0.0, <double*>C.data, &m)
    
    return C

方案2：利用zgemm间接实现（避免内存拷贝） #

如果可以接受使用复矩阵乘法接口zgemm而非dgemm，可以通过临时覆盖虚部的方式避免内存拷贝：

C语言示例（使用CBLAS） #

#include <complex.h>
#include <cblas.h>
#include <stdlib.h>

void multiply_complex_realparts_zgemm(int m, int k, int n, 
                                      double complex *A, double complex *B, double *C) {
    double *A_imag = malloc(m * k * sizeof(double));
    double *B_imag = malloc(k * n * sizeof(double));
    double complex *C_complex = malloc(m * n * sizeof(double complex));
    
    // 保存并清空A的虚部
    for (int j = 0; j < k; j++) {
        for (int i = 0; i < m; i++) {
            A_imag[i + j * m] = cimag(A[i + j * m]);
            A[i + j * m] = creal(A[i + j * m]) + 0.0*I;
        }
    }
    
    // 保存并清空B的虚部
    for (int j = 0; j < n; j++) {
        for (int t = 0; t < k; t++) {
            B_imag[t + j * k] = cimag(B[t + j * k]);
            B[t + j * k] = creal(B[t + j * k]) + 0.0*I;
        }
    }
    
    // 调用zgemm计算仅含实部的复矩阵乘法
    const char transa = 'N', transb = 'N';
    const double complex alpha = 1.0 + 0.0*I, beta = 0.0 + 0.0*I;
    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                m, n, k, &alpha, A, m, B, k, &beta, C_complex, m);
    
    // 提取结果的实部
    for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
            C[i + j * m] = creal(C_complex[i + j * m]);
        }
    }
    
    // 恢复原矩阵的虚部
    for (int j = 0; j < k; j++) {
        for (int i = 0; i < m; i++) {
            A[i + j * m] = creal(A[i + j * m]) + A_imag[i + j * m]*I;
        }
    }
    for (int j = 0; j < n; j++) {
        for (int t = 0; t < k; t++) {
            B[t + j * k] = creal(B[t + j * k]) + B_imag[t + j * k]*I;
        }
    }
    
    // 释放临时内存
    free(A_imag);
    free(B_imag);
    free(C_complex);
}

该方案无需拷贝实部，但需要临时修改原矩阵的虚部并恢复，实现复杂度更高，仅在内存极端紧张的场景下推荐使用。

内容的提问来源于stack exchange，提问作者G. Fougeron