Strassen's algorithm to compute matrix multiplication

//
//  main.cpp
//  Strassen
//
//  Created by Longxiang Lyu on 5/24/16.
//  Copyright (c) 2016 Longxiang Lyu. All rights reserved.
//

#include <iostream>
#include <vector>
#include <string>
#include <stdexcept>
#include <math.h>

using namespace std;

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

void zeroPadding(vector<vector<int>> &matrix)
{
    size_t sz = pow(2, (int)(sqrt(max(matrix.size(), matrix[0].size())) + 1));
    matrix.resize(sz);
    for (size_t i = 0; i < sz; ++i)
    {
        if (!matrix[i].empty())
            matrix[i].resize(sz);
        else
            matrix[i] = vector<int>(sz, 0);
    }
}

void sum(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &ret)
{
    size_t sz = A.size();

    for (int i = 0; i < sz; ++i)
        for (int j = 0; j < sz; ++j)
            ret[i][j] = (A[i][j] + B[i][j]);
}

void subtract(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &ret)
{
    size_t sz = A.size();
    // ret.clear();
    // ret.resize(sz);
    for (int i = 0; i < sz; ++i)
        for (int j = 0; j < sz; ++j)
            ret[i][j] = (A[i][j] - B[i][j]);
}



void strassenHelper(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &ret)
{
    if (A.size() == 1)
    {
        ret[0][0] = A[0][0] * B[0][0];
        return;
    }
    size_t sz = A.size();
    size_t new_sz = sz / 2;
    ret = vector<vector<int>>(sz, vector<int>(sz));

    vector<vector<int>> a11(new_sz), a12(new_sz), a21(new_sz), a22(new_sz), b11(new_sz), b12(new_sz), b21(new_sz), b22(new_sz);
    for (int i = 0; i < new_sz; ++i)
    {
        for (int j = 0; j < new_sz; ++j)
        {
            a11[i].push_back(A[i][j]);
            a12[i].push_back(A[i][j + new_sz]);
            a21[i].push_back(A[i + new_sz][j]);
            a22[i].push_back(A[i + new_sz][j + new_sz]);
            
            b11[i].push_back(B[i][j]);
            b12[i].push_back(B[i][j + new_sz]);
            b21[i].push_back(B[i + new_sz][j]);
            b22[i].push_back(B[i + new_sz][j + new_sz]);
        }
    }
    vector<vector<int>> result1(new_sz, vector<int>(new_sz, 0)), result2(new_sz, vector<int>(new_sz, 0));

    // p1
    vector<vector<int>> p1(new_sz, vector<int>(new_sz, 0));
    sum(a11, a22, result1);
    sum(b11, b22, result2);
    strassenHelper(result1, result2, p1);
    
    // p2
    vector<vector<int>> p2(new_sz, vector<int>(new_sz, 0));
    sum(a21, a22, result1);
    strassenHelper(result1, b11, p2);
    
    // p3
    vector<vector<int>> p3(new_sz, vector<int>(new_sz, 0));
    subtract(b12, b22, result2);
    strassenHelper(a11, result2, p3);
    
    // p4
    vector<vector<int>> p4(new_sz, vector<int>(new_sz, 0));
    subtract(b21, b11, result2);
    strassenHelper(a22, result2, p4);
    
    // p5
    vector<vector<int>> p5(new_sz, vector<int>(new_sz, 0));
    sum(a11, a12, result1);
    strassenHelper(result1, b22, p5);
    
    // p6
    vector<vector<int>> p6(new_sz, vector<int>(new_sz, 0));
    subtract(a21, a11, result1);
    sum(b11, b12, result2);
    strassenHelper(result1, result2, p6);
    
    // p7
    vector<vector<int>> p7(new_sz, vector<int>(new_sz, 0));
    subtract(a12, a22, result1);
    sum(b21, b22, result2);
    strassenHelper(result1, result2, p7);
    
    vector<vector<int>> c11(new_sz, vector<int>(new_sz, 0));
    vector<vector<int>> c12(new_sz, vector<int>(new_sz, 0));
    vector<vector<int>> c21(new_sz, vector<int>(new_sz, 0));
    vector<vector<int>> c22(new_sz, vector<int>(new_sz, 0));

    sum(p3, p5, c12);
    sum(p2, p4, c21);
    
    sum(p1, p4, result1);
    sum(result1, p7, result2);
    subtract(result2, p5, c11);
    
    sum(p1, p3, result1);
    sum(result1, p6, result2);
    subtract(result2, p2, c22);
    
    for (int i = 0; i < new_sz; ++i)
    {
        for (int j = 0; j < new_sz; ++j)
        {
            ret[i][j] = c11[i][j];
            ret[i][j + new_sz] = c12[i][j];
            ret[i + new_sz][j] = c21[i][j];
            ret[i + new_sz][j + new_sz] = c22[i][j];
        }
    }
    
}




void strassen(vector<vector<int>> &A, vector<vector<int>> &B, vector<vector<int>> &ret)
{
    if (A.empty() || B.empty())
        throw runtime_error("empty matrices");
    if (A[0].size() != B.size())
        throw runtime_error("A's col not equal B's row");
    zeroPadding(A);
    zeroPadding(B);
    strassenHelper(A, B, ret);
}

int main(int argc, const char * argv[]) {
    vector<vector<int>> A{{1, 2, 0}, {1, 2, 3}, {1, 2, 3}};
    vector<vector<int>> B{{1, 0, 1}, {1, 1, 1}, {2, 1, 1}};
    vector<vector<int>> ret(2, vector<int>(2));
    strassen(A, B, ret);
    printMatrix(ret);
    return 0;
}

Reference:

https://martin-thoma.com/strassen-algorithm-in-python-java-cpp/