Strassen快速矩阵乘法的C++实现

已于 2024-08-20 01:48:20 修改
阅读量367 收藏 6
点赞数 5
文章标签： c++ 矩阵 算法
于 2024-06-09 00:34:41 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/naturerun/article/details/139553395
版权

Strassen矩阵乘法是一个时间复杂度小于 O ( N 3 ) O(N^3) O(N3)的矩阵乘法算法，通过把朴素分治矩阵乘法的八次对子矩阵乘积计算的递归操作通过一系列优化技巧减少到七次，从而得到少于 O ( N 3 ) O(N^3) O(N3)的时间复杂度，该算法的详细讲解参见算法导论(第三版)中译本p43-p47这里只给出C++实现而不重复阐述
C++代码:

#include <iostream>
#include <vector>
#include <random>
using namespace std;

void matrixAdd(vector<vector<int>>& left, size_t left_up_x, size_t left_up_y, size_t left_down_x, size_t left_down_y,vector<vector<int>>& right, size_t right_up_x, size_t right_up_y, vector<vector<int>>& result)
{
    for (size_t i = 0; i <= left_down_y - left_up_y; ++i)
    {
        for (size_t j = 0; j <= left_down_x - left_up_x; ++j)
        {
            result[i][j] = left[i + left_up_y][j + left_up_x] + right[i + right_up_y][j + right_up_x];
        }
    }
}

void matrixSub(vector<vector<int>>& left, size_t left_up_x, size_t left_up_y, size_t left_down_x, size_t left_down_y, vector<vector<int>>& right, size_t right_up_x, size_t right_up_y, vector<vector<int>>& result)
{
    for (size_t i = 0; i <= left_down_y - left_up_y; ++i)
    {
        for (size_t j = 0; j <= left_down_x - left_up_x; ++j)
        {
            result[i][j] = left[i + left_up_y][j + left_up_x] - right[i + right_up_y][j + right_up_x];
        }
    }
}

void extract(vector<vector<int>>& left, size_t left_up_x, size_t left_up_y, size_t left_down_x, size_t left_down_y, vector<vector<int>>& result)
{
    for (size_t i = 0; i <= left_down_y - left_up_y; ++i)
    {
        for (size_t j = 0; j <= left_down_x - left_up_x; ++j)
        {
            result[i][j] = left[i + left_up_y][j + left_up_x];
        }
    }
}

void add(vector<vector<int>>& be_added, size_t left_up_x, size_t left_up_y, size_t left_down_x, size_t left_down_y, vector<vector<int>>& left)
{
    for (size_t i = 0; i + left_up_y <= left_down_y; ++i)
    {
        for (size_t j = 0; j + left_up_x <= left_down_x; ++j)
        {
            be_added[i + left_up_y][j + left_up_x] = left[i][j];
        }
    }
}

void Strassen(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &P)
{
    if (A.size() == 1 && B.size() == 1 && B[0].size() == 1)
    {
        P[0][0] = A[0][0] * B[0][0];
        return;
    }
    
    size_t s = A.size();
    size_t m = B[0].size();
    if (A.size() % 2 == 1)
    {
        A.push_back(vector<int>(B.size(), 0));
    }

    if (B.size() % 2 == 1)
    {
        for (size_t i = 0; i < A.size(); ++i)
        {
            A[i].push_back(0);
        }

        B.push_back(vector<int>(B[0].size(), 0));
    }

    if (B[0].size() % 2 == 1)
    {
        for (size_t j = 0; j < B.size(); ++j)
            B[j].push_back(0);
    }

    vector<vector<int>> S1(B.size()/2, vector<int>(B[0].size()/2, 0));
    matrixSub(B, B[0].size()/2, 0, B[0].size() - 1, B.size()/2 - 1, B, B[0].size()/2, B.size()/2, S1);
    vector<vector<int>> S2(A.size()/2, vector<int>(B.size()/2, 0));
    matrixAdd(A, 0, 0, B.size()/2 - 1, A.size()/2 - 1, A, B.size()/2, 0, S2);
    vector<vector<int>> S3(A.size() / 2, vector<int>(B.size() / 2, 0));
    matrixAdd(A, 0, A.size() / 2, B.size() / 2 - 1, A.size() - 1, A, B.size() / 2, A.size() / 2, S3);
    vector<vector<int>> S4(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixSub(B, 0, B.size()/2, B[0].size()/2 - 1, B.size() - 1, B, 0, 0, S4);
    vector<vector<int>> S5(A.size() / 2, vector<int>(B.size() / 2, 0));
    matrixAdd(A, 0, 0, B.size()/2 - 1, A.size()/2 - 1, A, B.size()/2, A.size()/2, S5);
    vector<vector<int>> S6(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(B, 0, 0, B[0].size()/2 - 1, B.size()/2 - 1, B, B[0].size()/2, B.size()/2, S6);
    vector<vector<int>> S7(A.size() / 2, vector<int>(B.size() / 2, 0));
    matrixSub(A, B.size()/2, 0, B.size()-1, A.size()/2 -1, A, B.size()/2, A.size()/2, S7);
    vector<vector<int>> S8(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(B, 0, B.size() / 2, B[0].size() / 2 - 1, B.size() - 1, B, B[0].size() / 2, B.size() / 2, S8);
    vector<vector<int>> S9(A.size() / 2, vector<int>(B.size() / 2, 0));
    matrixSub(A, 0, 0, B.size() / 2 - 1, A.size() / 2 - 1, A, 0, A.size() / 2, S9);
    vector<vector<int>> S10(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(B, 0, 0, B[0].size() / 2 - 1, B.size() / 2 - 1, B, B[0].size() / 2, 0, S10);

    vector<vector<int>> A11(A.size()/2, vector<int>(B.size()/2, 0));
    extract(A, 0, 0, B.size() / 2 - 1, A.size() / 2 - 1, A11);
    vector<vector<int>> P1(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(A11, S1, P1);

    vector<vector<int>> B22(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    extract(B, B[0].size()/2, B.size()/2, B[0].size() - 1, B.size() - 1, B22);
    vector<vector<int>> P2(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(S2, B22, P2);

    vector<vector<int>> B11(B.size() / 2, vector<int>(B[0].size() / 2, 0));
    extract(B, 0, 0, B[0].size() / 2 - 1, B.size() / 2 - 1, B11);
    vector<vector<int>> P3(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(S3, B11, P3);

    vector<vector<int>> A22(A.size() / 2, vector<int>(B.size() / 2, 0));
    extract(A, B.size() / 2, A.size() / 2, B.size() - 1, A.size() - 1, A22);
    vector<vector<int>> P4(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(A22, S4, P4);

    vector<vector<int>> P5(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(S5, S6, P5);

    vector<vector<int>> P6(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(S7, S8, P6);

    vector<vector<int>> P7(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    Strassen(S9, S10, P7);

    vector<vector<int>> C11(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(P5, 0, 0, P5[0].size() - 1, P5.size() - 1, P4, 0, 0, C11);
    matrixSub(C11, 0, 0, C11[0].size() - 1, C11.size() - 1, P2, 0, 0, C11);
    matrixAdd(C11, 0, 0, C11[0].size() - 1, C11.size() - 1, P6, 0, 0, C11);
    add(P, 0, 0, C11[0].size() - 1, C11.size() - 1, C11);

    vector<vector<int>> C12(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(P1, 0, 0, P1[0].size() - 1, P1.size() - 1, P2, 0, 0, C12);
    add(P, C11[0].size(), 0, P[0].size() - 1, C11.size() - 1, C12);

    vector<vector<int>> C21(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(P3, 0, 0, P3[0].size() - 1, P3.size() - 1, P4, 0, 0, C21);
    add(P, 0, C11.size(), C11[0].size() - 1, P.size() - 1, C21);

    vector<vector<int>> C22(A.size() / 2, vector<int>(B[0].size() / 2, 0));
    matrixAdd(P5, 0, 0, P5[0].size() - 1, P5.size() - 1, P1, 0, 0, C22);
    matrixSub(C22, 0, 0, C22[0].size() - 1, C22.size() - 1, P3, 0, 0, C22);
    matrixSub(C22, 0, 0, C22[0].size() - 1, C22.size() - 1, P7, 0, 0, C22);
    add(P, C11[0].size(), C11.size(), P[0].size() - 1, P.size() - 1, C22);
}

void matrixMultiply(vector<vector<int>>& left, vector<vector<int>>& right, vector<vector<int>>& result)
{
    for (size_t i = 0; i < left.size(); ++i)
    {
        for (size_t j = 0; j < right[0].size(); ++j)
        {
            for (size_t k = 0; k < right.size(); ++k)
                result[i][j] += left[i][k] * right[k][j];
        }
    }
}

int main()
{
    for (int NA = 1; NA <= 10; ++NA)
    {
        for (int MB = 1; MB <= 10; ++MB)
        {
            for (int MA = 1; MA <= 10; ++MA)
            {
                int NB = MA;
                vector<int> testA(NA * MA, 0);
                for (int i = 0; i < NA * MA; ++i)
                    testA[i] = i + 1;
                shuffle(testA.begin(), testA.end(), default_random_engine());
                vector<vector<int>> A(NA, vector<int>(MA, 0));

                size_t k = 0;
                for (size_t i = 0; i < A.size(); ++i)
                {
                    for (size_t j = 0; j < A[i].size(); ++j)
                    {
                        A[i][j] = testA[k++];
                    }
                }

                vector<int> testB(NB * MB, 0);
                for (int i = 0; i < NB * MB; ++i)
                    testB[i] = i + 1;
                shuffle(testB.begin(), testB.end(), default_random_engine());
                vector<vector<int>> B(NB, vector<int>(MB, 0));

                k = 0;
                for (size_t i = 0; i < B.size(); ++i)
                {
                    for (size_t j = 0; j < B[i].size(); ++j)
                    {
                        B[i][j] = testB[k++];
                    }
                }

                cout << "相乘的两个矩阵为" << endl;
                for (size_t i = 0; i < A.size(); ++i)
                {
                    for (size_t j = 0; j < A[i].size(); ++j)
                    {
                        cout << A[i][j] << " ";
                    }
                    cout << endl;
                }
                cout << endl;

                for (size_t i = 0; i < B.size(); ++i)
                {
                    for (size_t j = 0; j < B[i].size(); ++j)
                    {
                        cout << B[i][j] << " ";
                    }
                    cout << endl;
                }

                vector<vector<int>> result1(NA, vector<int>(MB, 0));
                vector<vector<int>> result2(NA, vector<int>(MB, 0));
                vector<vector<int>> _B = B;
                vector<vector<int>> _A = A;
                Strassen(A, B, result1);
                matrixMultiply(_A, _B, result2);
                for (size_t i = 0; i < result1.size(); ++i)
                {
                    for (size_t j = 0; j < result1[i].size(); ++j)
                    {
                        if (result1[i][j] != result2[i][j])
                        {
                            cout << "矩阵相乘结果错误!" << endl;
                            exit(-1);
                        }
                    }
                }
                cout << "矩阵相乘结果正确,乘积为" << endl;
                for (size_t i = 0; i < result1.size(); ++i)
                {
                    for (size_t j = 0; j < result1[i].size(); ++j)
                    {
                        cout << result1[i][j] << " ";
                    }
                    cout << endl;
                }
            }
        }
    }

   /* const int NA = 3;
    const int MA = 4;
    const int NB = 4;
    const int MB = 5;
    vector<int> testA(NA * MA, 0);
    for (int i = 0; i < NA * MA; ++i)
        testA[i] = i + 1;
    shuffle(testA.begin(), testA.end(), default_random_engine());
    vector<vector<int>> A(NA, vector<int>(MA, 0));

    size_t k = 0;
    for (size_t i = 0; i < A.size(); ++i)
    {
        for (size_t j = 0; j < A[i].size(); ++j)
        {
            A[i][j] = testA[k++];
        }
    }

    vector<int> testB(NB * MB, 0);
    for (int i = 0; i < NB * MB; ++i)
        testB[i] = i + 1;
    shuffle(testB.begin(), testB.end(), default_random_engine());
    vector<vector<int>> B(NB, vector<int>(MB, 0));

    k = 0;
    for (size_t i = 0; i < B.size(); ++i)
    {
        for (size_t j = 0; j < B[i].size(); ++j)
        {
            B[i][j] = testB[k++];
        }
    }

    cout << "相乘的两个矩阵为" << endl;
    for (size_t i = 0; i < A.size(); ++i)
    {
        for (size_t j = 0; j < A[i].size(); ++j)
        {
            cout << A[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;

    for (size_t i = 0; i < B.size(); ++i)
    {
        for (size_t j = 0; j < B[i].size(); ++j)
        {
            cout << B[i][j] << " ";
        }
        cout << endl;
    }

    vector<vector<int>> result1(NA, vector<int>(MB, 0));
    vector<vector<int>> result2(NA, vector<int>(MB, 0));
    vector<vector<int>> _B = B;
    vector<vector<int>> _A = A;
    Strassen(A, B, result1);
    matrixMultiply(_A, _B, result2);
    for (size_t i = 0; i < result1.size(); ++i)
    {
        for (size_t j = 0; j < result1[i].size(); ++j)
        {
            if (result1[i][j] != result2[i][j])
            {
                cout << "矩阵相乘结果错误!" << endl;
                exit(-1);
            }
        }
    }
    cout << "矩阵相乘结果正确,乘积为" << endl;
    for (size_t i = 0; i < result1.size(); ++i)
    {
        for (size_t j = 0; j < result1[i].size(); ++j)
        {
            cout << result1[i][j] << " ";
        }
        cout << endl;
    }*/
    return 0;
}
naturerun
关注
Strassen矩阵乘法 c++实现
04-11
算法分析与设计课程作业，单个cpp文件
Strassen矩阵乘法—分治法实现代码
12-13
Strassen矩阵乘法是一种高效的矩阵乘法算法，由德国数学家Volker Strassen在1969年提出，是分治策略的一个经典应用。相比于传统的矩阵乘法算法，它的主要优势在于减少了运算次数，尤其在处理大尺寸矩阵时，效率提升...
C语言strassen矩阵乘法代码
04-05
strassen矩阵乘法的C代码 【问题描述】 从文件arr.in中读入一个m行k列的整数矩阵a和一个k行n列的整数矩阵b(1 , k, n )，在标准输出上输出这两个矩阵的乘积。 【输入形式】 输入文件arr.in中有m+k行，前m...
矩阵乘法的strassen算法(C++实现)
04-23
一般情况下矩阵乘法需要三个for循环，时间复杂度为O(n^3)，现在我们将矩阵分块如图：( 来自MIT算法导论 ) 一般算法需要八次乘法 r = a * e + b * g ; s = a * f + b * h ; t = c * e + d * g; u = c * f + d * h; ...
Strassen矩阵乘法
11-28
Strassen矩阵乘法是一种高效的矩阵乘法算法，由德国数学家Volker Strassen在1969年提出，是计算理论中的一个重要成果。相比于传统的矩阵乘法算法，它通过分治策略来降低运算次数，尤其是在处理大矩阵时，效率优势...
juzhen.zip_矩阵乘法 快速
09-24
在实际编程实现中，如压缩包中的"juzhen.cpp"文件，可能会包含一个矩阵乘法的C++实现，可能采用了Strassen算法或者其他优化过的分治策略。代码通常包括矩阵的定义、分解、乘法和合并操作，以及递归调用。优化的重点...
C++ Strassen算法代码的实现
08-19
Strassen算法的基本思想是将矩阵分割成小块，然后使用递归的方式来实现矩阵乘法。通过这种方式，Strassen算法可以将矩阵乘法的时间复杂度降低到O(n^2.81)。 在实现Strassen算法时，我们需要定义一个矩阵结构体来...
Desktop_C++_算法导论_矩阵乘法_递归_
10-03
在本项目中，我们主要探讨的是C++编程语言中的矩阵乘法算法，特别是基于《算法导论》第六版第四章77页所描述的递归方法。递归是一种强大的编程技术，它允许函数调用自身来解决问题。在这个特定的场景中，递归被用来...
TestMatrix_多种C++矩阵乘法算法_avx_
09-29
在本文中，我们将深入探讨标题为“TestMatrix_多种C++矩阵乘法算法_avx_”的项目，该项目旨在实现并比较多种C++矩阵乘法算法，特别关注使用Intel AVX（Advanced Vector Extensions）指令集的优化技术。矩阵乘法是...
Strassen:施特拉森矩阵乘法
07-11
总结来说，Strassen矩阵乘法是算法设计中的一个经典例子，它展示了如何通过分治策略和创新性思考来改进传统算法。C++实现Strassen算法可以帮助我们更好地理解和掌握这种高效计算技术，并为解决其他复杂问题提供思路...
快速矩阵乘法_C++_MATLAB_下载.zip
04-13
本资源包"快速矩阵乘法_C++_MATLAB_下载.zip"提供了C++和MATLAB两种编程语言实现快速矩阵乘法的方法。接下来，我们将深入探讨矩阵乘法的基本概念、快速矩阵乘法算法以及C++和MATLAB的实现细节。 首先，矩阵乘法是...
Strassen 算法 C++ 实现 任意矩阵相乘
11-18
Strassen 算法 C++ 实现 任意矩阵相乘 用Command line 输入编制好的两个矩阵，输出相乘的结果矩阵。如果需要手动输入，删除程序里的相关语句，添加输入命令即可。
Strassen矩阵 分治算法 递归 C++
07-28
Strassen是采用分治算法的思想，将所给矩阵分成2阶矩阵 分治的方法循序渐进处理各个小矩阵的相乘，一个矩阵可以分成更多小的矩阵的。
Strassen矩阵乘法Karatsuba快速乘法最近点对问题源代码(C/C++程序)
05-19
在提供的压缩文件中，可能包含C或C++语言编写的程序，用于实现Strassen矩阵乘法和Karatsuba快速乘法。这些程序可能包含递归函数，用于执行分治过程，并展示了如何在实际编程环境中应用这些理论算法。理解并运行这些...
Strassen矩阵相乘算法的c++代码实现.docx
11-25
在本篇文章中，我们将深入探讨Strassen矩阵乘法的原理及其C++实现细节。 #### 二、Strassen算法的基本思想 1. **递归分割**：将大矩阵分解为若干个小矩阵。 2. **计算子矩阵**：通过特定方式计算7个较小规模的矩阵...
蓝桥杯国赛题之C++矩阵乘法.zip
最新发布
06-03
在C++中实现矩阵乘法，我们可以使用二维数组来表示矩阵，并采用双重循环来计算每个元素。以下是一个基本的矩阵乘法函数实现： ```cpp #include using namespace std; const int MAX_SIZE = 100; // 假设矩阵的...
aaa.zip_矩阵乘法
09-19
在实现上，有多种方法可以执行矩阵乘法，包括串行乘法、并行乘法（如OpenMP、CUDA）、Strassen算法、Coppersmith-Winograd算法等。串行乘法是最直观的实现方式，但效率相对较低，尤其是在处理大规模矩阵时。512*512...
juzhengchengfa.rar_矩阵乘法
09-24
**C++实现矩阵乘法** 在C++中，我们可以使用二维数组来表示矩阵，然后通过两层循环进行乘法运算。以下是实现的基本步骤： 1. **定义矩阵**：首先，定义两个二维数组代表矩阵A和B，以及一个二维数组C用于存储乘积...
分治法实现矩阵相乘
03-13
这个过程可以通过递归算法实现，每次将问题规模减半，直到矩阵的大小为1x1，此时矩阵相乘变为直接的乘法运算。然后逐级回溯，将子问题的解合并得到最终结果。分治法的递归公式如下： \[ C = \begin{bmatrix} C_{11}...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值