使用 gsl 的复杂矩阵乘法

Question

我有一个非常简单的问题——我只想在 C 语言的 gsl 中将两个具有复杂条目的矩阵相乘。例如，我想定义一个函数

gsl_matrix_complex *multiply( gsl_matrix_complex *A, gsl_matrix_complex *B ) {
???
}

我已经看到当条目是双倍时完成的，但我无法找出适当的概括。任何帮助将不胜感激!

再次感谢。

Answer 1

要在 GSL 中乘以矩阵，您需要使用稍微神秘的 BLAS interface .

假设您要相乘的矩阵没有任何已知结构(例如对称或厄密等)。要使用的适当函数是 gsl_blas_zgemm。函数名中看似随意的一串字母的翻译是:

z  = double precision complex
ge = general matrices
mm = matrix-matrix multiplication

(有关所有缩写的列表，请参阅 here)

使用文档中的信息，我们可以将您的函数编写为

#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>

gsl_matrix_complex *multiply( gsl_matrix_complex *A, gsl_matrix_complex *B)
{
    gsl_matrix_complex *result = gsl_matrix_complex_alloc(A->size1, B->size2);

    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans,
                   GSL_COMPLEX_ONE, A, B,
                   GSL_COMPLEX_ZERO, result);

    return result;
}

请注意，与其在乘法函数中分配结果矩阵，不如传入一个已在调用函数中分配的结果矩阵，例如使用签名:

void multiply*(gsl_matrix_complex *A, gsl_matrix_complex *B, gsl_matrix_complex *result);

然后省略分配步骤和显式返回。这样做的原因是它允许您最大限度地减少代码执行的堆分配次数，并使您更难意外编写泄漏内存的代码。