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;
}