#include <iostream>
|
|
#include <cstdlib>
|
|
#include <random>
|
|
#include <vector>
|
|
#include <boost/thread.hpp>
|
|
#include <boost/bind.hpp>
|
|
#include <boost/chrono.hpp>
|
|
/*g++ -std=c++11 fox_boost.cpp -o fox_boost -lboost_thread -lboost_chrono -lboost_system
|
|
./fox_boost 8 2
|
|
Note that this implementation:
|
|
|
|
Uses shared memory (all threads access the full matrices)
|
|
|
|
Requires careful synchronization to avoid race conditions
|
|
|
|
May not be as efficient as MPI for large matrices due to shared memory contention
|
|
*/
|
|
|
|
using namespace std;
|
|
namespace chrono = boost::chrono;
|
|
|
|
// Matrix structure to hold data and dimensions
|
|
struct Matrix {
|
|
int n;
|
|
double** data;
|
|
|
|
Matrix(int size) : n(size) {
|
|
data = new double*[n];
|
|
for (int i = 0; i < n; ++i) {
|
|
data[i] = new double[n](); // Initialize to zero
|
|
}
|
|
}
|
|
|
|
~Matrix() {
|
|
for (int i = 0; i < n; ++i) {
|
|
delete[] data[i];
|
|
}
|
|
delete[] data;
|
|
}
|
|
};
|
|
|
|
// Initialize matrix with random values
|
|
void init_matrix(Matrix& mat, bool zero = false) {
|
|
random_device rd;
|
|
mt19937 gen(rd());
|
|
uniform_real_distribution<double> dis(0.0, 1.0);
|
|
|
|
for (int i = 0; i < mat.n; ++i) {
|
|
for (int j = 0; j < mat.n; ++j) {
|
|
mat.data[i][j] = zero ? 0.0 : dis(gen);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Print matrix (for debugging)
|
|
void print_matrix(const Matrix& mat, const string& name) {
|
|
cout << name << ":\n";
|
|
for (int i = 0; i < mat.n; ++i) {
|
|
for (int j = 0; j < mat.n; ++j) {
|
|
cout << mat.data[i][j] << " ";
|
|
}
|
|
cout << endl;
|
|
}
|
|
cout << endl;
|
|
}
|
|
|
|
// Worker thread for Fox algorithm
|
|
class FoxWorker {
|
|
int block_size;
|
|
int grid_size;
|
|
int row, col;
|
|
Matrix& a;
|
|
Matrix& b;
|
|
Matrix& c;
|
|
boost::barrier& sync_barrier;
|
|
vector<boost::barrier*>& row_barriers;
|
|
|
|
public:
|
|
FoxWorker(int bs, int gs, int r, int c, Matrix& ma, Matrix& mb, Matrix& mc,
|
|
boost::barrier& sb, vector<boost::barrier*>& rb)
|
|
: block_size(bs), grid_size(gs), row(r), col(c),
|
|
a(ma), b(mb), c(mc), sync_barrier(sb), row_barriers(rb) {}
|
|
|
|
void operator()() {
|
|
// Calculate block boundaries
|
|
int row_start = row * block_size;
|
|
int row_end = row_start + block_size;
|
|
int col_start = col * block_size;
|
|
int col_end = col_start + block_size;
|
|
|
|
// Temporary storage for broadcast block
|
|
Matrix temp_block(block_size);
|
|
|
|
for (int step = 0; step < grid_size; ++step) {
|
|
// Determine broadcast root for this step
|
|
int bcast_root = (row + step) % grid_size;
|
|
|
|
// Broadcast phase - only one thread per row participates
|
|
if (col == bcast_root) {
|
|
// Copy my block of A to temp_block
|
|
for (int i = 0; i < block_size; ++i) {
|
|
for (int j = 0; j < block_size; ++j) {
|
|
temp_block.data[i][j] = a.data[row_start + i][col_start + j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Wait for all threads in this row to reach broadcast point
|
|
row_barriers[row]->wait();
|
|
|
|
// All threads in row now have access to temp_block
|
|
|
|
// Local multiplication phase
|
|
for (int i = 0; i < block_size; ++i) {
|
|
for (int j = 0; j < block_size; ++j) {
|
|
for (int k = 0; k < block_size; ++k) {
|
|
c.data[row_start + i][col_start + j] +=
|
|
temp_block.data[i][k] * b.data[row_start + k][col_start + j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Shift blocks of B upwards
|
|
if (row > 0) {
|
|
// Send my block of B to the neighbor above
|
|
for (int i = 0; i < block_size; ++i) {
|
|
for (int j = 0; j < block_size; ++j) {
|
|
double temp = b.data[row_start + i][col_start + j];
|
|
b.data[row_start + i][col_start + j] = b.data[(row_start + i - block_size + a.n) % a.n][col_start + j];
|
|
b.data[(row_start + i - block_size + a.n) % a.n][col_start + j] = temp;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Synchronize all threads before next step
|
|
sync_barrier.wait();
|
|
}
|
|
}
|
|
};
|
|
|
|
// Fox algorithm using Boost Threads
|
|
void fox_multiply(Matrix& a, Matrix& b, Matrix& c, int grid_size) {
|
|
int block_size = a.n / grid_size;
|
|
|
|
// Create synchronization barriers
|
|
boost::barrier sync_barrier(grid_size * grid_size);
|
|
vector<boost::barrier*> row_barriers(grid_size);
|
|
for (int i = 0; i < grid_size; ++i) {
|
|
row_barriers[i] = new boost::barrier(grid_size);
|
|
}
|
|
|
|
// Create worker threads
|
|
boost::thread_group workers;
|
|
for (int i = 0; i < grid_size; ++i) {
|
|
for (int j = 0; j < grid_size; ++j) {
|
|
workers.create_thread(FoxWorker(block_size, grid_size, i, j,
|
|
a, b, c, sync_barrier, row_barriers));
|
|
}
|
|
}
|
|
|
|
// Wait for all threads to complete
|
|
workers.join_all();
|
|
|
|
// Clean up barriers
|
|
for (int i = 0; i < grid_size; ++i) {
|
|
delete row_barriers[i];
|
|
}
|
|
}
|
|
|
|
int main(int argc, char** argv) {
|
|
if (argc != 3) {
|
|
cerr << "Usage: " << argv[0] << " <matrix_size> <grid_size>" << endl;
|
|
return 1;
|
|
}
|
|
|
|
int n = atoi(argv[1]);
|
|
int grid_size = atoi(argv[2]);
|
|
if (n <= 0 || grid_size <= 0) {
|
|
cerr << "Matrix size and grid size must be positive" << endl;
|
|
return 1;
|
|
}
|
|
if (n % grid_size != 0) {
|
|
cerr << "Matrix size must be divisible by grid size" << endl;
|
|
return 1;
|
|
}
|
|
|
|
chrono::system_clock::time_point setup_start, setup_end;
|
|
chrono::system_clock::time_point comp_start, comp_end;
|
|
|
|
setup_start = chrono::system_clock::now();
|
|
|
|
Matrix a(n), b(n), c(n);
|
|
init_matrix(a);
|
|
init_matrix(b);
|
|
init_matrix(c, true); // Initialize to zero
|
|
|
|
setup_end = chrono::system_clock::now();
|
|
|
|
if (n <= 8) {
|
|
print_matrix(a, "Matrix A");
|
|
print_matrix(b, "Matrix B");
|
|
}
|
|
|
|
comp_start = chrono::system_clock::now();
|
|
|
|
fox_multiply(a, b, c, grid_size);
|
|
|
|
comp_end = chrono::system_clock::now();
|
|
|
|
if (n <= 8) {
|
|
print_matrix(c, "Result Matrix C");
|
|
}
|
|
|
|
auto setup_time = chrono::duration_cast<chrono::milliseconds>(setup_end - setup_start);
|
|
auto comp_time = chrono::duration_cast<chrono::milliseconds>(comp_end - comp_start);
|
|
|
|
cout << "n,grid_size,setup(ms),exec(ms)" << endl;
|
|
cout << n << "," << grid_size << "," << setup_time.count() << "," << comp_time.count() << endl;
|
|
|
|
return 0;
|
|
}
|