#include <immintrin.h>

#include <string.h>
#include "sgemm_intrinsics.h"

///Compute the Matrix-Matrix product of the x,y matrices
/**	@param[out] matOut : result
 * 	@param matX : left matrix
 * 	@param matY : right matrix
 * 	@param size : size of the square matrices
*/
void sgemm(float* matOut, const float* matX, const float* matY, long unsigned int size){
	memset(matOut, 0, sizeof(float)*size*size);
	long unsigned int vecSize(VECTOR_ALIGNEMENT/sizeof(float));
	long unsigned int nbVec(size/vecSize);
	
	for(long unsigned int i(0lu); i < size; ++i){			//Iterate over X rows
		const float * rowX = matX + i*size;				//Get current X row
		float * rowOut = matOut + i*size;					//Get current out row
		for(long unsigned int k(0lu); k < size; ++k){					//Part of dot product
			__m256 regX = _mm256_broadcast_ss(rowX + k);
			
			const float* rowY = matY + k*size;				//Get current Y row
			for(long unsigned int j(0lu); j < nbVec; ++j){
				__m256 regY = _mm256_load_ps(rowY + vecSize*j);
				
				__m256 regRes = _mm256_load_ps(rowOut + vecSize*j);
				
				regRes = _mm256_add_ps(regRes, _mm256_mul_ps(regX, regY));
				
				_mm256_store_ps(rowOut + vecSize*j, regRes);
			}
		}
	}
}