矩阵相乘strassen-c++代码实现及运行实例结果

理论知识参见算法导论第三版以及百度


伪代码



c++代码

#include <iostream>

using namespace std;

#define N 2//以二维方阵为例

template<typename T>//使用模板保证矩阵可以为不同类型如int,double等
void output(T D[N][N],int n);

template<typename T>
void strassen(T A[N][N],T B[N][N],T C[N][N],int n);

int main()
{
    int A[N][N]={1,2,3,4};
    cout<<"矩阵A的值"<<endl;
    output(A,N);//先输出矩阵A的值
    int B[N][N]={5,6,7,8};
    cout<<"矩阵B的值"<<endl;
    output(B,N);//先输出矩阵B的值
    int C[N][N];
    strassen(A,B,C,N);
    cout<<"矩阵C的值"<<endl;
    output(C,N);//输出矩阵C的值

    return 0;
}

template<typename T>
void output(T D[N][N],int n)
{
    for(int i=0;i<n;++i)
    {
        for(int j=0;j<n;++j)
        {
            cout<<D[i][j]<<" ";
        }
        cout<<endl;
    }


}

template<typename T>
void matrixAdd(T a[N][N],T b[N][N],T c[N][N],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];
}

template<typename T>
void matrixSub(T a[N][N],T b[N][N],T c[N][N],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];
}

template<typename T>
void matrixMul(T a[N][N],T b[N][N],T c[N][N])//定义矩阵乘法
{
    for(int i=0;i<2;++i)//小于2是因为strassen采用递归,递归结束标志是最终分成二阶矩阵
        for(int j=0;j<2;++j)
        {
            c[i][j]=0;
            for(int k=0;k<2;++k)
                c[i][j]+=a[i][k]*b[k][j];
        }
}

template<typename T>
void strassen(T A[N][N],T B[N][N],T C[N][N],int n)
{
    T A11[N][N],A12[N][N],A21[N][N],A22[N][N];//将矩阵分块
    T B11[N][N],B12[N][N],B21[N][N],B22[N][N];
    T C11[N][N],C12[N][N],C21[N][N],C22[N][N];
    T S1[N][N],S2[N][N],S3[N][N],S4[N][N],S5[N][N],S6[N][N],S7[N][N];//由strassen定义的7个系数
    T temp1[N][N],temp2[N][N];//存储中间量
    if(n==2)//递归结束标志
        matrixMul(A,B,C);
    else
    {
        for(int i=0;i<n/2;++i)
            for(int j=0;j<n/2;++j)
            {
                A11[i][j]=A[i][j];//将矩阵分成相同的4个块
                A12[i][j]=A[i][j+n/2];
                A21[i][j]=A[i+n/2][j];
                A22[i][j]=A[i+n/2][j+n/2];

                B11[i][j]=B[i][j];
                B12[i][j]=B[i][j+n/2];
                B21[i][j]=B[i+n/2][j];
                B22[i][j]=B[i+n/2][j+n/2];
            }
     /*S1=(A11+A22)×(B11+B22)*/
     matrixAdd(A11,A22,temp1,n/2);
     matrixAdd(B11,B22,temp2,n/2);
     strassen(temp1,temp2,S1,n/2);
     /*S2=(A21+A22)×B11*/
     matrixAdd(A21,A22,temp1,n/2);
     strassen(temp1,B11,S2,n/2);
     /*S3=A11*(B12+B22)*/
     matrixSub(B12,B22,temp1,n/2);
     strassen(A11,temp1,S3,n/2);
     /*S4=A22×(B21+B11)*/
     matrixSub(B21,B11,temp1,n/2);
     strassen(A22,temp1,S4,n/2);
     /*S5=(A11+A12)×B22*/
     matrixAdd(A11,A12,temp1,n/2);
     strassen(temp1,B22,S5,n/2);
     /*S6=(A11+A22)×(B11+B22)*/
     matrixSub(A21,A11,temp1,n/2);
     matrixAdd(B11,B12,temp2,n/2);
     strassen(temp1,temp2,S6,n/2);
     /*S7=(A12+A22)×(B21+B22)*/
     matrixSub(A12,A22,temp1,n/2);
     matrixAdd(B21,B22,temp2,n/2);
     strassen(temp1,temp2,S7,n/2);
     /*C11 = S1+S4-S5+S7*/
     matrixAdd(S1,S4,temp1,n/2);
     matrixSub(S7,S5,temp2,n/2);
     matrixAdd(temp1,temp2,C11,n/2);
     /*C12 = S3+S5*/
     matrixAdd(S3,S5,C12,n/2);
     /*C21 = S2+S4*/
     matrixAdd(S2,S4,C21,n/2);
     /*C22 = S1-S2+S3+S6*/
     matrixSub(S1,S2,temp1,n/2);
     matrixAdd(S3,S6,temp2,n/2);
     matrixAdd(temp1,temp2,C22,n/2);

     for(int i=0;i<n/2;++i)//逆转化
         for(int j=0;j<n/2;++j)
         {
             C[i][j]=C11[i][j];
             C[i][j+n/2]=C12[i][j];
             C[i+n/2][j]=C21[i][j];
             C[i+n/2][j+n/2]= C22[i][j];
         }
    }
}

运行结果




  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
以下是一个简单的 Strassen 矩阵乘法的分治法的 C 代码实现: ```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 时,使用普通矩阵乘法。分割矩阵使用了动态内存分配,临时矩阵的内存同样使用了动态分配,并在最后释放了所有动态分配的内存。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值