SIMD Matrix

SimdMatrix.cpp

#include <xmmintrin.h>
#include <iostream>
#include <iomanip>
#include <random>
#include <vector>
#include <sstream>
#include <chrono>
using namespace std;

union alignas(16) Float16
{
	float values[16];
	__m128 groups[4];
};

ostream& WriteRowMajor(ostream& stream, const Float16& rmi)
{
	for (int i = 0; i < 4; ++i)
	{
		for (int j = 0; j < 4; ++j)
		{
			stream << setw(4) << rmi.values[i * 4 + j];
		}

		stream << '\n';
	}

	return stream;
}

ostream& WriteColumnMajor(ostream& stream, const Float16& cmi)
{
	for (int i = 0; i < 4; ++i)
	{
		for (int j = 0; j < 4; ++j)
		{
			stream << setw(4) << cmi.values[j * 4 + i];
		}

		stream << '\n';
	}

	return stream;
}

void Multiply(const Float16& rmi, const Float16& cmi, Float16& rmo)
{
	for (int i = 0; i < 4; ++i)
	{
		for (int j = 0; j < 4; ++j)
		{
			float total = 0.0f;

			for (int k = 0; k < 4; ++k)
				total += rmi.values[i * 4 + k] * cmi.values[j * 4 + k];

			rmo.values[i * 4 + j] = total;
		}
	}
}

void SimdMultiply(const Float16& rmi, const Float16& cmi, Float16& rmo)
{
	for (int i = 0; i < 4; ++i)
	{
		for (int j = 0; j < 4; ++j)
		{
			__m128 products = _mm_mul_ps(rmi.groups[i], cmi.groups[j]);

			float total = 0.0f;

			for (int k = 0; k < 4; ++k)
				total += ((float*)&products)[k];

			rmo.values[i * 4 + j] = total;
		}
	}
}

void Randomize(Float16& matrix, mt19937_64& mt)
{
	uniform_int_distribution<int> distribution(-9, 9);
	for (int i = 0; i < 16; ++i) matrix.values[i] = distribution(mt);
}

int main(int argc, char** argv)
{
	int matrixCount = 4;

	if (argc > 1)
	{
		stringstream ss;
		ss << argv[1];
		ss >> matrixCount;
	}

	mt19937_64 mt;
	vector<Float16> matrices(matrixCount * 4);
	Float16* rowMajorInputs = matrices.data();
	Float16* columnMajorInputs = rowMajorInputs + matrixCount;
	Float16* rowMajorOutputs = columnMajorInputs + matrixCount;
	Float16* rowMajorOutputsSimd = rowMajorOutputs + matrixCount;

	cout << "Generating data..." << endl;
	for (int i = 0; i < matrixCount; ++i)
	{
		Randomize(rowMajorInputs[i], mt);
		Randomize(columnMajorInputs[i], mt);
	}

	cout << "Simple matrix multiplication...";
	cout.flush();
	auto start = chrono::system_clock::now();
	for (int i = 0; i < matrixCount; ++i)
		Multiply(rowMajorInputs[i], columnMajorInputs[i], rowMajorOutputs[i]);
	auto end = chrono::system_clock::now();
	chrono::duration<double> elapsed = end - start;
	cout << elapsed.count() << endl;

	cout << "SIMD matrix multiplication...";
	cout.flush();
	start = chrono::system_clock::now();
	for (int i = 0; i < matrixCount; ++i)
		SimdMultiply(rowMajorInputs[i], columnMajorInputs[i], rowMajorOutputsSimd[i]);	
	end = chrono::system_clock::now();
	elapsed = end - start;
	cout << elapsed.count() << endl;

	WriteRowMajor(cout, rowMajorOutputs[0]) << '\n';
	WriteRowMajor(cout, rowMajorOutputsSimd[0]);

	return 0;
}