You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
182 lines
5.2 KiB
182 lines
5.2 KiB
/*Strassen 1D - Matrix sizes at the power of to to subdivice properly*/
|
|
/*All matrices use 1D arrays with flattened indexing.
|
|
|
|
Recursion is threaded for the 7 Strassen multiplications up to threshold 64 using boost::bind*/
|
|
#include <iostream>
|
|
#include <cstdlib>
|
|
#include <cstring>
|
|
#include <boost/thread.hpp>
|
|
#include <boost/bind/bind.hpp>
|
|
#include <boost/chrono.hpp>
|
|
#include <random>
|
|
using namespace std;
|
|
using namespace boost::placeholders;
|
|
|
|
double* create_matrix(int size) {
|
|
return new double[size * size]();
|
|
}
|
|
|
|
void cleanup_matrix(double* matrix) {
|
|
delete[] matrix;
|
|
}
|
|
|
|
void add(double* a, double* b, double* result, int size) {
|
|
for (int i = 0; i < size * size; ++i)
|
|
result[i] = a[i] + b[i];
|
|
}
|
|
|
|
void subtract(double* a, double* b, double* result, int size) {
|
|
for (int i = 0; i < size * size; ++i)
|
|
result[i] = a[i] - b[i];
|
|
}
|
|
|
|
void naive_mul(double* a, double* b, double* c, int size) {
|
|
for (int i = 0; i < size; ++i)
|
|
for (int j = 0; j < size; ++j)
|
|
for (int k = 0; k < size; ++k)
|
|
c[i * size + j] += a[i * size + k] * b[k * size + j];
|
|
}
|
|
|
|
void strassen(double* a, double* b, double* c, int size, int threshold = 64) {
|
|
if (size <= threshold) {
|
|
naive_mul(a, b, c, size);
|
|
return;
|
|
}
|
|
|
|
int newSize = size / 2;
|
|
int block = newSize * newSize;
|
|
|
|
double *a11 = a;
|
|
double *a12 = a + newSize;
|
|
double *a21 = a + newSize * size;
|
|
double *a22 = a + newSize * size + newSize;
|
|
|
|
double *b11 = b;
|
|
double *b12 = b + newSize;
|
|
double *b21 = b + newSize * size;
|
|
double *b22 = b + newSize * size + newSize;
|
|
|
|
double *c11 = c;
|
|
double *c12 = c + newSize;
|
|
double *c21 = c + newSize * size;
|
|
double *c22 = c + newSize * size + newSize;
|
|
|
|
// Allocate temp matrices
|
|
double *m1 = create_matrix(newSize), *m2 = create_matrix(newSize), *m3 = create_matrix(newSize);
|
|
double *m4 = create_matrix(newSize), *m5 = create_matrix(newSize), *m6 = create_matrix(newSize), *m7 = create_matrix(newSize);
|
|
double *aTemp = create_matrix(newSize), *bTemp = create_matrix(newSize);
|
|
|
|
boost::thread_group group;
|
|
|
|
// M1 = (A11 + A22) * (B11 + B22)
|
|
add(a11, a22, aTemp, newSize);
|
|
add(b11, b22, bTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, aTemp, bTemp, m1, newSize, threshold));
|
|
|
|
// M2 = (A21 + A22) * B11
|
|
add(a21, a22, aTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, aTemp, b11, m2, newSize, threshold));
|
|
|
|
// M3 = A11 * (B12 - B22)
|
|
subtract(b12, b22, bTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, a11, bTemp, m3, newSize, threshold));
|
|
|
|
// M4 = A22 * (B21 - B11)
|
|
subtract(b21, b11, bTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, a22, bTemp, m4, newSize, threshold));
|
|
|
|
// M5 = (A11 + A12) * B22
|
|
add(a11, a12, aTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, aTemp, b22, m5, newSize, threshold));
|
|
|
|
// M6 = (A21 - A11) * (B11 + B12)
|
|
subtract(a21, a11, aTemp, newSize);
|
|
add(b11, b12, bTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, aTemp, bTemp, m6, newSize, threshold));
|
|
|
|
// M7 = (A12 - A22) * (B21 + B22)
|
|
subtract(a12, a22, aTemp, newSize);
|
|
add(b21, b22, bTemp, newSize);
|
|
group.create_thread(boost::bind(strassen, aTemp, bTemp, m7, newSize, threshold));
|
|
|
|
group.join_all();
|
|
|
|
// C11 = M1 + M4 - M5 + M7
|
|
for (int i = 0; i < block; ++i)
|
|
c11[i] = m1[i] + m4[i] - m5[i] + m7[i];
|
|
|
|
// C12 = M3 + M5
|
|
for (int i = 0; i < block; ++i)
|
|
c12[i] = m3[i] + m5[i];
|
|
|
|
// C21 = M2 + M4
|
|
for (int i = 0; i < block; ++i)
|
|
c21[i] = m2[i] + m4[i];
|
|
|
|
// C22 = M1 - M2 + M3 + M6
|
|
for (int i = 0; i < block; ++i)
|
|
c22[i] = m1[i] - m2[i] + m3[i] + m6[i];
|
|
|
|
cleanup_matrix(m1);
|
|
cleanup_matrix(m2);
|
|
cleanup_matrix(m3);
|
|
cleanup_matrix(m4);
|
|
cleanup_matrix(m5);
|
|
cleanup_matrix(m6);
|
|
cleanup_matrix(m7);
|
|
cleanup_matrix(aTemp);
|
|
cleanup_matrix(bTemp);
|
|
}
|
|
|
|
void init_matrix(double* mat, int n, bool zero = false) {
|
|
random_device rd;
|
|
mt19937 gen(rd());
|
|
uniform_real_distribution<double> dis(0.0, 1.0);
|
|
for (int i = 0; i < n * n; ++i)
|
|
mat[i] = zero ? 0.0 : dis(gen);
|
|
}
|
|
|
|
int main(int argc, char** argv) {
|
|
if (argc != 2) {
|
|
cerr << "Usage: " << argv[0] << " <matrix_size (power of 2)>" << endl;
|
|
return 1;
|
|
}
|
|
|
|
int n = atoi(argv[1]);
|
|
if ((n & (n - 1)) != 0) {
|
|
cerr << "Matrix size must be a power of 2!" << endl;
|
|
return 1;
|
|
}
|
|
|
|
boost::chrono::system_clock::time_point t1, t2, t3, t4;
|
|
|
|
t1 = boost::chrono::system_clock::now();
|
|
|
|
double* a = create_matrix(n);
|
|
double* b = create_matrix(n);
|
|
double* c = create_matrix(n);
|
|
|
|
init_matrix(a, n);
|
|
init_matrix(b, n);
|
|
init_matrix(c, n, true);
|
|
|
|
t2 = boost::chrono::system_clock::now();
|
|
t3 = boost::chrono::system_clock::now();
|
|
|
|
strassen(a, b, c, n);
|
|
|
|
t4 = boost::chrono::system_clock::now();
|
|
|
|
auto setup_time = boost::chrono::duration_cast<boost::chrono::milliseconds>(t2 - t1);
|
|
auto exec_time = boost::chrono::duration_cast<boost::chrono::milliseconds>(t4 - t3);
|
|
|
|
cout << "n,setup(ms),exec(ms)" << endl;
|
|
cout << n << "," << setup_time.count() << "," << exec_time.count() << endl;
|
|
|
|
cleanup_matrix(a);
|
|
cleanup_matrix(b);
|
|
cleanup_matrix(c);
|
|
|
|
return 0;
|
|
}
|
|
|