如何使用此 C 代码使用 Strassen 算法将两个矩阵相乘？

Question

我一直在寻找 C 中 Strassen's Algorithm 的实现，我在最后找到了这段代码。

使用multiply函数：

void multiply(int n, matrix a, matrix b, matrix c, matrix d);

将两个矩阵a、b 相乘并将结果放入c（d 是一个中间矩阵）。矩阵a 和b 应具有以下类型：

typedef union _matrix 
{
    double **d;
    union _matrix **p;
} *matrix;

我已经动态分配了四个矩阵a、b、c、d（二维数组）并将它们的地址分配给字段_matrix.d：

#include "strassen.h"

#define SIZE 50 

int main(int argc, char *argv[])
{
    double ** matA, ** matB, ** matC, ** matD;
    union _matrix ma, mb, mc, md; 
    int i = 0, j = 0, n;

    matA = (double **) malloc(sizeof(double *) * SIZE);
    for (i = 0; i < SIZE; i++)
        matA[i] = (double *) malloc(sizeof(double) * SIZE); 
    // Do the same for matB, matC, matD.

    ma.d = matA;
    mb.d = matB;
    mc.d = matC;
    md.d = matD;

    // Initialize matC and matD to 0.

    // Read n.

    // Read matA and matB.

    multiply(n, &ma, &mb, &mc, &md);
    return 0;
}

此代码成功编译，但由于n > BREAK 而崩溃。

strassen.c：

#include "strassen.h"

/* c = a * b */
void multiply(int n, matrix a, matrix b, matrix c, matrix d)
{
    if (n <= BREAK) {
      double sum, **p = a->d, **q = b->d, **r = c->d;
      int i, j, k;

      for (i = 0; i < n; i++)
         for (j = 0; j < n; j++) {
            for (sum = 0., k = 0; k < n; k++)
               sum += p[i][k] * q[k][j];
            r[i][j] = sum;
         }
    } else {
        n /= 2;
        sub(n, a12, a22, d11);
        add(n, b21, b22, d12);
        multiply(n, d11, d12, c11, d21);
        sub(n, a21, a11, d11);
        add(n, b11, b12, d12);
        multiply(n, d11, d12, c22, d21);
        add(n, a11, a12, d11);
        multiply(n, d11, b22, c12, d12);
        sub(n, c11, c12, c11);
        sub(n, b21, b11, d11);
        multiply(n, a22, d11, c21, d12);
        add(n, c21, c11, c11);
        sub(n, b12, b22, d11);
        multiply(n, a11, d11, d12, d21);
        add(n, d12, c12, c12);
        add(n, d12, c22, c22);
        add(n, a21, a22, d11);
        multiply(n, d11, b11, d12, d21);
        add(n, d12, c21, c21);
        sub(n, c22, d12, c22);
        add(n, a11, a22, d11);
        add(n, b11, b22, d12);
        multiply(n, d11, d12, d21, d22);
        add(n, d21, c11, c11);
        add(n, d21, c22, c22);
    }
}

/* c = a + b */
void add(int n, matrix a, matrix b, matrix c)
{
    if (n <= BREAK) {
        double **p = a->d, **q = b->d, **r = c->d;
        int i, j;

        for (i = 0; i < n; i++)
           for (j = 0; j < n; j++)
              r[i][j] = p[i][j] + q[i][j];
    } else {
        n /= 2;
        add(n, a11, b11, c11);
        add(n, a12, b12, c12);
        add(n, a21, b21, c21);
        add(n, a22, b22, c22);
    }
}

/* c = a - b */
void sub(int n, matrix a, matrix b, matrix c)
{
    if (n <= BREAK) {
        double **p = a->d, **q = b->d, **r = c->d;
        int i, j;

        for (i = 0; i < n; i++)
           for (j = 0; j < n; j++)
              r[i][j] = p[i][j] - q[i][j];
    } else {
        n /= 2;
        sub(n, a11, b11, c11);
        sub(n, a12, b12, c12);
        sub(n, a21, b21, c21);
        sub(n, a22, b22, c22);
    }
}

strassen.h：

#define BREAK 8   

typedef union _matrix {
    double **d;
    union _matrix **p;
} *matrix;

/* Notational shorthand to access submatrices for matrices named a, b, c, d */

#define a11 a->p[0]
#define a12 a->p[1]
#define a21 a->p[2]
#define a22 a->p[3]
#define b11 b->p[0]
#define b12 b->p[1]
#define b21 b->p[2]
#define b22 b->p[3]
#define c11 c->p[0]
#define c12 c->p[1]
#define c21 c->p[2]
#define c22 c->p[3]
#define d11 d->p[0]
#define d12 d->p[1]
#define d21 d->p[2]
#define d22 d->p[3]

我的问题是如何使用函数multiply（如何实现矩阵）。

strassen.h

strassen.c

Answer 1

与Atom said 一样，您需要为两个矩阵正确初始化matrix.p。

1) 首先，您的matrix 是union，所以p 本质上变成d 解释为_matrix **，这在这里没有意义——这就是它崩溃的原因。您可能需要将matrix 改为struct。
最后，p 根据定义是一个子矩阵数组，因此它应该是一个struct _matrix *（并且您需要在需要时使用malloc 实际数组）或struct _matrix[4]（这是不可能的:)）。

typedef struct _matrix 
{
    double **d;
    struct _matrix *p;
} *matrix;

2) 现在，让我们看看p 应该是什么。

                           │
A.d ->  d1 -> a[1,1] a[1,2]│a[1,3] a[1,4]
        d2 -> a[2,1] a[2,2]│a[2,3] a[2,4]
             ─────────────────────────────
        d3 -> a[3,1] a[3,2]│a[3,3] a[3,4]
        d4 -> a[4,1] a[4,2]│a[4,3] a[4,4]
                           │

p 指向matrix 结构的数组！特殊之处在于使这些结构中的d 指向A内部 这样(p[k].d)[i][j] 就是相应子矩阵的元素：

p[0].d -> p01 -> a[1,1]    p[1].d -> p11 -> a[1,3]
          p02 -> a[2,1]              p12 -> a[2,3]

p[2].d -> p21 -> a[3,1]    p[3].d -> p31 -> a[3,3]
          p22 -> a[4,1]              p32 -> a[4,3]

你现在能推导出算法来为任意偶数大小的正方形A初始化p吗？

首先什么时候初始化它？ ;)

Answer 2

当 n > BREAK 时，矩阵乘法算法使用分层矩阵表示（union _matrix 的字段 p，而不是字段 d）。

在分配内存和初始化矩阵a 和b 时，您需要针对分层表示调整代码。