strassen矩阵乘法分治法C++代码实现

以下是一个简单的 Strassen 矩阵乘法的分治法的 C 代码实现：

#include <stdio.h>
#include <stdlib.h>

// 普通矩阵乘法
void matmul(int **a, int **b, int **c, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            c[i][j] = 0;
            for (int k = 0; k < n; k++) {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
}

// 矩阵加法
void matadd(int **a, int **b, int **c, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            c[i][j] = a[i][j] + b[i][j];
        }
    }
}

// 矩阵减法
void matsub(int **a, int **b, int **c, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            c[i][j] = a[i][j] - b[i][j];
        }
    }
}

// Strassen 矩阵乘法
void strassen(int **a, int **b, int **c, int n) {
    // 当矩阵大小小于等于 64 时，使用普通矩阵乘法
    if (n <= 64) {
        matmul(a, b, c, n);
        return;
    }

    // 分割矩阵
    int m = n / 2;
    int **a11 = (int **) malloc(m * sizeof(int *));
    int **a12 = (int **) malloc(m * sizeof(int *));
    int **a21 = (int **) malloc(m * sizeof(int *));
    int **a22 = (int **) malloc(m * sizeof(int *));
    int **b11 = (int **) malloc(m * sizeof(int *));
    int **b12 = (int **) malloc(m * sizeof(int *));
    int **b21 = (int **) malloc(m * sizeof(int *));
    int **b22 = (int **) malloc(m * sizeof(int *));
    int **c11 = (int **) malloc(m * sizeof(int *));
    int **c12 = (int **) malloc(m * sizeof(int *));
    int **c21 = (int **) malloc(m * sizeof(int *));
    int **c22 = (int **) malloc(m * sizeof(int *));
    for (int i = 0; i < m; i++) {
        a11[i] = (int *) malloc(m * sizeof(int));
        a12[i] = (int *) malloc(m * sizeof(int));
        a21[i] = (int *) malloc(m * sizeof(int));
        a22[i] = (int *) malloc(m * sizeof(int));
        b11[i] = (int *) malloc(m * sizeof(int));
        b12[i] = (int *) malloc(m * sizeof(int));
        b21[i] = (int *) malloc(m * sizeof(int));
        b22[i] = (int *) malloc(m * sizeof(int));
        c11[i] = (int *) malloc(m * sizeof(int));
        c12[i] = (int *) malloc(m * sizeof(int));
        c21[i] = (int *) malloc(m * sizeof(int));
        c22[i] = (int *) malloc(m * sizeof(int));
    }
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            a11[i][j] = a[i][j];
            a12[i][j] = a[i][j + m];
            a21[i][j] = a[i + m][j];
            a22[i][j] = a[i + m][j + m];
            b11[i][j] = b[i][j];
            b12[i][j] = b[i][j + m];
            b21[i][j] = b[i + m][j];
            b22[i][j] = b[i + m][j + m];
        }
    }

    // 临时矩阵
    int **p1 = (int **) malloc(m * sizeof(int *));
    int **p2 = (int **) malloc(m * sizeof(int *));
    int **p3 = (int **) malloc(m * sizeof(int *));
    int **p4 = (int **) malloc(m * sizeof(int *));
    int **p5 = (int **) malloc(m * sizeof(int *));
    int **p6 = (int **) malloc(m * sizeof(int *));
    int **p7 = (int **) malloc(m * sizeof(int *));
    for (int i = 0; i < m; i++) {
        p1[i] = (int *) malloc(m * sizeof(int));
        p2[i] = (int *) malloc(m * sizeof(int));
        p3[i] = (int *) malloc(m * sizeof(int));
        p4[i] = (int *) malloc(m * sizeof(int));
        p5[i] = (int *) malloc(m * sizeof(int));
        p6[i] = (int *) malloc(m * sizeof(int));
        p7[i] = (int *) malloc(m * sizeof(int));
    }

    // 计算七个子问题
    matsub(b12, b22, p1, m); // p1 = b12 - b22
    strassen(a11, p1, p2, m); // p2 = a11 x p1
    matadd(a11, a12, p1, m); // p1 = a11 + a12
    strassen(p1, b22, p3, m); // p3 = p1 x b22
    matadd(a21, a22, p1, m); // p1 = a21 + a22
    strassen(p1, b11, p4, m); // p4 = p1 x b11
    matsub(b21, b11, p1, m); // p1 = b21 - b11
    strassen(a22, p1, p5, m); // p5 = a22 x p1
    matadd(a11, a22, p1, m); // p1 = a11 + a22
    matadd(b11, b22, p2, m); // p2 = b11 + b22
    strassen(p1, p2, p6, m); // p6 = p1 x p2
    matadd(p4, p5, p1, m); // p1 = p4 + p5
    matadd(p1, p6, p2, m); // p2 = p1 + p6
    matsub(p2, p3, c11, m); // c11 = p2 - p3
    matadd(p4, p1, c12, m); // c12 = p4 + p1
    matadd(p5, p3, c21, m); // c21 = p5 + p3
    matsub(p6, p2, c22, m); // c22 = p6 - p2

    // 合并结果
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            c[i][j] = c11[i][j];
            c[i][j + m] = c12[i][j];
            c[i + m][j] = c21[i][j];
            c[i + m][j + m] = c22[i][j];
        }
    }

    // 释放内存
    for (int i = 0; i < m; i++) {
        free(a11[i]);
        free(a12[i]);
        free(a21[i]);
        free(a22[i]);
        free(b11[i]);
        free(b12[i]);
        free(b21[i]);
        free(b22[i]);
        free(c11[i]);
        free(c12[i]);
        free(c21[i]);
        free(c22[i]);
        free(p1[i]);
        free(p2[i]);
        free(p3[i]);
        free(p4[i]);
        free(p5[i]);
        free(p6[i]);
        free(p7[i]);
    }
    free(a11);
    free(a12);
    free(a21);
    free(a22);
    free(b11);
    free(b12);
    free(b21);
    free(b22);
    free(c11);
    free(c12);
    free(c21);
    free(c22);
    free(p1);
    free(p2);
    free(p3);
    free(p4);
    free(p5);
    free(p6);
    free(p7);
}

int main() {
    int n = 8;
    int **a = (int **) malloc(n * sizeof(int *));
    int **b = (int **) malloc(n * sizeof(int *));
    int **c = (int **) malloc(n * sizeof(int *));
    for (int i = 0; i < n; i++) {
        a[i] = (int *) malloc(n * sizeof(int));
        b[i] = (int *) malloc(n * sizeof(int));
        c[i] = (int *) malloc(n * sizeof(int));
    }
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            a[i][j] = rand() % 10;
            b[i][j] = rand() % 10;
            c[i][j] = 0;
        }
    }

    strassen(a, b, c, n);

    printf("Matrix A:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%d ", a[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    printf("Matrix B:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%d ", b[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    printf("Matrix C:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%d ", c[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    for (int i = 0; i < n; i++) {
        free(a[i]);
        free(b[i]);
        free(c[i]);
    }
    free(a);
    free(b);
    free(c);

    return 0;
}

此代码中，实现了普通矩阵乘法、矩阵加法、矩阵减法和 Strassen 矩阵乘法。在 Strassen 矩阵乘法中，当矩阵大小小于等于 64 时，使用普通矩阵乘法。分割矩阵使用了动态内存分配，临时矩阵的内存同样使用了动态分配，并在最后释放了所有动态分配的内存。