Stranssen算法的实现与优化

通过本实验，比较并验证了普通矩阵算法和Stranssen算法的运行效率及时间复杂度，并对Stranssen算法进行了改进，使之适应非2^n的情况
 

#include <iostream>
#include <vector>
#include <ctime>
#include <cstdlib>
#include <chrono>

using namespace std;
using namespace std::chrono;

vector<vector<int>> generateRandomMatrix(int N, int s, int t) {
    vector<vector<int>> matrix(N, vector<int>(N));
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            matrix[i][j] = rand() % (t - s + 1) + s;
        }
    }
    return matrix;
}

void printMatrix(const vector<vector<int>>& matrix) {
    for (const auto& row : matrix) {
        for (int val : row) {
            cout << val << " ";
        }
        cout << endl;
    }
}

vector<vector<int>> subtractMatrices(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> result(n, vector<int>(n));

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            result[i][j] = A[i][j] - B[i][j];
        }
    }

    return result;
}

vector<vector<int>> addMatrices(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> result(n, vector<int>(n));

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            result[i][j] = A[i][j] + B[i][j];
        }
    }

    return result;
}

vector<vector<int>> matrixMultiplication(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> result(n, vector<int>(n, 0));

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            for (int k = 0; k < n; ++k) {
                result[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return result;
}

vector<vector<int>> strassenMatrixMultiplication(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> result(n, vector<int>(n, 0));

    if (n == 1) {
        result[0][0] = A[0][0] * B[0][0];
    } else {
        int half = n / 2;

        vector<vector<int>> A11(half, vector<int>(half));
        vector<vector<int>> A12(half, vector<int>(half));
        vector<vector<int>> A21(half, vector<int>(half));
        vector<vector<int>> A22(half, vector<int>(half));
        vector<vector<int>> B11(half, vector<int>(half));
        vector<vector<int>> B12(half, vector<int>(half));
        vector<vector<int>> B21(half, vector<int>(half));
        vector<vector<int>> B22(half, vector<int>(half));

        // Divide A and B into 4 sub-matrices
        for (int i = 0; i < half; ++i) {
            for (int j = 0; j < half; ++j) {
                A11[i][j] = A[i][j];
                A12[i][j] = A[i][j + half];
                A21[i][j] = A[i + half][j];
                A22[i][j] = A[i + half][j + half];

                B11[i][j] = B[i][j];
                B12[i][j] = B[i][j + half];
                B21[i][j] = B[i + half][j];
                B22[i][j] = B[i + half][j + half];
            }
        }

        vector<vector<int>> P1 = strassenMatrixMultiplication(A11, subtractMatrices(B12, B22));
        vector<vector<int>> P2 = strassenMatrixMultiplication(addMatrices(A11, A12), B22);
        vector<vector<int>> P3 = strassenMatrixMultiplication(addMatrices(A21, A22), B11);
        vector<vector<int>> P4 = strassenMatrixMultiplication(A22, subtractMatrices(B21, B11));
        vector<vector<int>> P5 = strassenMatrixMultiplication(addMatrices(A11, A22), addMatrices(B11, B22));
        vector<vector<int>> P6 = strassenMatrixMultiplication(subtractMatrices(A12, A22), addMatrices(B21, B22));
        vector<vector<int>> P7 = strassenMatrixMultiplication(subtractMatrices(A11, A21), addMatrices(B11, B12));

        vector<vector<int>> C11 = addMatrices(subtractMatrices(addMatrices(P5, P4), P2), P6);
        vector<vector<int>> C12 = addMatrices(P1, P2);
        vector<vector<int>> C21 = addMatrices(P3, P4);
        vector<vector<int>> C22 = subtractMatrices(subtractMatrices(addMatrices(P5, P1), P3), P7);

        for (int i = 0; i < half; ++i) {
            for (int j = 0; j < half; ++j) {
                result[i][j] = C11[i][j];
                result[i][j + half] = C12[i][j];
                result[i + half][j] = C21[i][j];
                result[i + half][j + half] = C22[i][j];
            }
        }
    }

    return result;
}

int main() {
    srand(time(0));

    
    vector<int> sizes = {1024, 2048,4096,8192,16384}; 
    int s = 1; 
    int t = 10; 

    for (int N : sizes) {
        
        vector<vector<int>> A = generateRandomMatrix(N, s, t);
        vector<vector<int>> B = generateRandomMatrix(N, s, t);

       
        auto start = high_resolution_clock::now();
        matrixMultiplication(A, B);
        auto stop = high_resolution_clock::now();
        auto duration = duration_cast<microseconds>(stop - start);
        cout << "Normal Matrix Multiplication for size " << N << ": " << duration.count() << " microseconds" << endl;

        
        start = high_resolution_clock::now();
        strassenMatrixMultiplication(A, B);
        stop = high_resolution_clock::now();
        duration = duration_cast<microseconds>(stop - start);
        cout << "Strassen's Matrix Multiplication for size " << N << ": " << duration.count() << " microseconds" << endl;
    }

    return 0;
}

改进版
再main函数前加入补0的函数

void pad_matrix(int n, int **original, int **padded) {
    int padded_n = 1;
    while (padded_n < n) {
        padded_n *= 2;    
}
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            padded[i][j] = original[i][j];
        }
    }
    for (int i = n; i < padded_n; i++) {
        for (int j = 0; j < padded_n; j++) {
            padded[i][j] = 0;
        }
    }
    for (int i = 0; i < padded_n; i++) {
        for (int j = n; j < padded_n; j++) {
            padded[i][j] = 0;
        }
    }
}