c++マトリクス乗算とブロックマトリクス乗算を実現

マトリクスAサイズ:m*p、マトリクスBサイズ:p*m、結果マトリクスCサイズ:m*n、ブロックサイズk*k.
くだらないことは言わないで、原理も言わないで、直接コードに行きます

#include "iostream"
#include "random"
#include "vector"

using namespace std;

vector> initMatrix(int m, int n);

vector> multiply(vector> &matrix_a, vector> &matrix_b);

vector>
multiply_block(vector> &matrix_a, vector> &matrix_b, int block_size);


bool isEqual(vector> &matrix_a, vector> &matrix_b);

int main() {
    //m * p, p * n, and k represent the size of block
    int m = 300, p = 300, n = 300, k = 7;
    auto matrix_a = initMatrix(m, p), matrix_b = initMatrix(p, n);
    vector> matrix_c(m, vector(n, 0));
    auto matrix_c_1 = multiply_block(matrix_a, matrix_b, k);
    auto matrix_c_2 = multiply(matrix_a, matrix_b);
    cout << isEqual(matrix_c_1, matrix_c_2);
    return 0;
}

vector> initMatrix(int m, int n) {
    vector> res(m, vector(n, 0));
    default_random_engine engine;
    uniform_real_distribution u(0.0, 200);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            res[i][j] = u(engine);
        }
    }
    return res;
}


// normal matrix multiplication 
vector> multiply(vector> &matrix_a, vector> &matrix_b) {
    int m = matrix_a.size(), p = matrix_a[0].size(), n = matrix_b[0].size();
    vector> res(m, vector(n, 0));
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            for (int k = 0; k < p; ++k) {
                res[i][j] += matrix_a[i][k] * matrix_b[k][j];
            }
        }
    }
    return res;
}

// block-based matrix multiplication
vector>
multiply_block(vector> &matrix_a, vector> &matrix_b, int block_size) {
    int m = matrix_a.size(), p = matrix_a[0].size(), n = matrix_b[0].size();
    vector> res(m, vector(n, 0));
    for (int bi = 0; bi < m; bi += block_size) {
        for (int bj = 0; bj < n; bj += block_size) {
            for (int bk = 0; bk < p; bk += block_size) {
                for (int i = bi; i < min(bi + block_size, m); ++i) {
                    for (int j = bj; j < min(bj + block_size, n); ++j) {
                        for (int k = bk; k < min(bk + block_size, p); ++k) {
                            res[i][j] += matrix_a[i][k] * matrix_b[k][j];
                        }
                    }
                }
            }
        }
    }
    return res;
}

// to check whether the two results are same
bool isEqual(vector> &matrix_a, vector> &matrix_b) {
    int m_a = matrix_a.size(), n_a = matrix_a[0].size();
    int m_b = matrix_b.size(), n_b = matrix_b[0].size();
    if (m_a != m_b || n_a != n_b)
        return false;
    for (int i = 0; i < m_a; ++i) {
        for (int j = 0; j < n_a; ++j) {
            if (matrix_a[i][j] != matrix_b[i][j])
                return false;
        }
    }
    return true;
}