Source codes

27 files changed, 5017 insertions, 0 deletions

diff --git a/Whisper/CPU/BufferAllocator.cpp b/Whisper/CPU/BufferAllocator.cpp
new file mode 100644
index 0000000..156382c
--- /dev/null
+++ b/Whisper/CPU/BufferAllocator.cpp
@@ -0,0 +1,145 @@
+#include <stdafx.h>
+#include "BufferAllocator.h"
+#include <immintrin.h>
+#include <ammintrin.h>
+using namespace CpuCompute;
+
+HRESULT BufferAllocator::create( size_t cb )
+{
+	CHECK( buffer.allocate( cb ) );
+	head = 0;
+	size = cb;
+	dbgMarkUninitializedMemory( buffer.pointer(), cb );
+	return S_OK;
+}
+
+namespace
+{
+	// Round up the integer by 32 bytes
+	__forceinline size_t roundUpAlloc( size_t cb )
+	{
+		const size_t mask = 31;
+		cb += mask;
+		// We require AVX1+FMA3 support, might as well use BMI1
+		return _andn_u64( mask, cb );
+	}
+}
+
+void* BufferAllocator::allocate( size_t cb, size_t align ) noexcept
+{
+	assert( align <= 32 );
+	cb = roundUpAlloc( cb );
+
+	uint8_t* pointer = buffer.pointer();
+	if( head + cb > size || nullptr == pointer )
+	{
+		logError( u8"BufferAllocator.allocate, not enough capacity" );
+		return nullptr;
+	}
+
+	void* const res = pointer + head;
+	head += cb;
+	assert( head <= size );
+	dbgMarkUninitializedMemory( res, cb );
+	return res;
+}
+
+namespace
+{
+	// 2 MB of memory, we hope the OS kernel will then be smart enough to give us large pages.
+	constexpr size_t virtualAllocGranularityExp2 = 21;
+
+	constexpr size_t virtualAllocGranularityMask = ( ( (size_t)1 ) << virtualAllocGranularityExp2 ) - 1;
+
+	// Round up the integer by 2 megabytes
+	__forceinline size_t roundUpVirtualAlloc( size_t cb )
+	{
+		const size_t mask = virtualAllocGranularityMask;
+		cb += mask;
+		return _andn_u64( mask, cb );
+	}
+}
+
+HRESULT VirtualAllocator::create( size_t cb )
+{
+	if( nullptr != pointer )
+		return HRESULT_FROM_WIN32( ERROR_ALREADY_INITIALIZED );
+	cb = roundUpVirtualAlloc( cb );
+	pointer = (uint8_t*)VirtualAlloc( NULL, cb, MEM_RESERVE, PAGE_READWRITE );
+	if( nullptr != pointer )
+	{
+		head = 0;
+		sizeAllocated = 0;
+		sizeVirtual = cb;
+		return S_OK;
+	}
+
+	const HRESULT hr = getLastHr();
+	logErrorHr( hr, u8"VirtualAlloc failed" );
+	return hr;
+}
+
+void* VirtualAllocator::allocate( size_t cb, size_t align ) noexcept
+{
+	assert( align <= 32 );
+	cb = roundUpAlloc( cb );
+
+	const size_t newHead = head + cb;
+	if( newHead <= sizeAllocated )
+	{
+		void* const res = pointer + head;
+		head = newHead;
+		dbgMarkUninitializedMemory( res, cb );
+		return res;
+	}
+
+	if( newHead <= sizeVirtual )
+	{
+		uint8_t* const ptrCommit = pointer + sizeAllocated;
+		const size_t cbCommit = roundUpVirtualAlloc( newHead ) - sizeAllocated;
+		void* const res = VirtualAlloc( ptrCommit, cbCommit, MEM_COMMIT, PAGE_READWRITE );
+		if( nullptr != res )
+		{
+			sizeAllocated += cbCommit;
+			assert( sizeAllocated <= sizeVirtual );
+			void* const res = pointer + head;
+			head = newHead;
+			dbgMarkUninitializedMemory( res, cb );
+			return res;
+		}
+
+		const HRESULT hr = getLastHr();
+		logErrorHr( hr, u8"VirtualAllocator.allocate, VirtualAlloc failed" );
+		return nullptr;
+	}
+
+	logError( u8"VirtualAllocator.allocate, not enough arena capacity" );
+	return nullptr;
+}
+
+VirtualAllocator::~VirtualAllocator()
+{
+	if( nullptr == pointer )
+		return;
+
+	if( VirtualFree( pointer, 0, MEM_RELEASE ) )
+	{
+		pointer = nullptr;
+		return;
+	}
+
+	const HRESULT hr = getLastHr();
+	logErrorHr( hr, u8"VirtualFree failed" );
+}
+
+#ifndef NDEBUG
+// Reusing Microsoft's magic numbers: https://asawicki.info/news_1292_magic_numbers_in_visual_c
+void CpuCompute::dbgMarkUninitializedMemory( void* pv, size_t cb )
+{
+	__stosb( (uint8_t*)pv, 0xCD, cb );
+}
+void CpuCompute::dbgMarkFreedMemory( void* pv, size_t cb )
+{
+	__stosd( (DWORD*)pv, 0xFEEEFEEEu, cb / 4 );
+}
+#endif
\ No newline at end of file
diff --git a/Whisper/CPU/BufferAllocator.h b/Whisper/CPU/BufferAllocator.h
new file mode 100644
index 0000000..750a565
--- /dev/null
+++ b/Whisper/CPU/BufferAllocator.h
@@ -0,0 +1,64 @@
+#pragma once
+#include "LargeBuffer.h"
+#include "Tensor.h"
+
+namespace CpuCompute
+{
+#ifdef NDEBUG
+	inline void dbgMarkUninitializedMemory( void* pv, size_t cb ) { }
+	inline void dbgMarkFreedMemory( void* pv, size_t cb ) { }
+#else
+	void dbgMarkUninitializedMemory( void* pv, size_t cb );
+	void dbgMarkFreedMemory( void* pv, size_t cb );
+#endif
+
+	// An implementation of arena allocator which slices pieces of a large buffer allocated in advance
+	class BufferAllocator : public iArenaAllocator
+	{
+		LargeBuffer buffer;
+		size_t head = 0;
+		size_t size = 0;
+
+		void resetArena() noexcept override final
+		{
+			head = 0;
+			dbgMarkFreedMemory( buffer.pointer(), size );
+		}
+
+		void* allocate( size_t cb, size_t align ) noexcept override final;
+
+	public:
+		BufferAllocator() = default;
+		BufferAllocator( const BufferAllocator& ) = delete;
+		~BufferAllocator() = default;
+
+		// Allocate a large buffer with the specified count of bytes
+		HRESULT create( size_t cb );
+	};
+
+	// An implementation of arena allocator which allocates a large chunk of virtual memory, and maps new physical pages into that memory region as needed.
+	class VirtualAllocator : public iArenaAllocator
+	{
+		uint8_t* pointer = nullptr;
+		size_t head = 0;
+		size_t sizeAllocated = 0;
+		size_t sizeVirtual = 0;
+
+		void resetArena() noexcept override final
+		{
+			head = 0;
+			dbgMarkFreedMemory( pointer, sizeAllocated );
+		}
+
+		void* allocate( size_t cb, size_t align ) noexcept override final;
+
+	public:
+
+		VirtualAllocator() = default;
+		VirtualAllocator( const VirtualAllocator& ) = delete;
+		~VirtualAllocator();
+
+		// Reserve virtual memory space for the specified count of bytes in the arena, but don't allocate any pages
+		HRESULT create( size_t cb );
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/DecoderTensors.cpp b/Whisper/CPU/DecoderTensors.cpp
new file mode 100644
index 0000000..22de476
--- /dev/null
+++ b/Whisper/CPU/DecoderTensors.cpp
@@ -0,0 +1,68 @@
+#include "stdafx.h"
+#include "DecoderTensors.h"
+using namespace CpuCompute;
+
+#if TENSOR_GGML_COMPAT
+namespace
+{
+	class CompatContext
+	{
+		std::vector<ggml_tensor>& vec;
+		size_t index;
+
+	public:
+		CompatContext( std::vector<ggml_tensor>& dest, size_t layers ) :
+			vec( dest )
+		{
+			constexpr size_t tensorsPerLayer = 21;
+			const size_t count = tensorsPerLayer * layers + 4;
+			vec.resize( count );
+			index = 0;
+		}
+
+		void add( const Tensor& rsi, ggml_tensor*& res )
+		{
+			ggml_tensor& ten = vec[ index ];
+			index++;
+			ten = rsi.ggml();
+			res = &ten;
+		}
+
+		void add2( const TensorPair& rsi, ggml_tensor*& w, ggml_tensor*& b )
+		{
+			add( rsi.w, w );
+			add( rsi.b, b );
+		}
+
+		bool isComplete() const
+		{
+			return index == vec.size();
+		}
+	};
+}
+
+void DecoderTensors::makeCompatTensors()
+{
+	CompatContext ctx( ggml, layers.size() );
+
+	ctx.add( positionalEmbedding, d_pe );
+	ctx.add( tokenEmbedding, d_te );
+	ctx.add2( ln, d_ln_w, d_ln_b );
+
+	for( auto& i : layers )
+	{
+		ctx.add2( i.attnLn0, i.attn_ln_0_w, i.attn_ln_0_b );
+		ctx.add2( i.attnLn1, i.attn_ln_1_w, i.attn_ln_1_b );
+		ctx.add2( i.attnQuery, i.attn_q_w, i.attn_q_b );
+		ctx.add( i.attnKey, i.attn_k_w );
+		ctx.add2( i.attnValue, i.attn_v_w, i.attn_v_b );
+		ctx.add2( i.crossAttnLn0, i.cross_attn_ln_0_w, i.cross_attn_ln_0_b );
+		ctx.add2( i.crossAttnLn1, i.cross_attn_ln_1_w, i.cross_attn_ln_1_b );
+		ctx.add2( i.crossAttnQuery, i.cross_attn_q_w, i.cross_attn_q_b );
+		ctx.add2( i.mlpLn, i.mlp_ln_w, i.mlp_ln_b );
+		ctx.add2( i.mlp0, i.mlp_0_w, i.mlp_0_b );
+		ctx.add2( i.mlp1, i.mlp_1_w, i.mlp_1_b );
+	}
+	assert( ctx.isComplete() );
+}
+#endif
\ No newline at end of file
diff --git a/Whisper/CPU/DecoderTensors.h b/Whisper/CPU/DecoderTensors.h
new file mode 100644
index 0000000..2efa519
--- /dev/null
+++ b/Whisper/CPU/DecoderTensors.h
@@ -0,0 +1,131 @@
+#pragma once
+#include <vector>
+#include "Tensor.h"
+#include "LargeBuffer.h"
+#if TENSOR_GGML_COMPAT
+#include "../source/ggml.h"
+#endif
+
+namespace CpuCompute
+{
+	// A set of tensors for one decoder's layer
+	struct LayerDecoder
+	{
+		// decoder.blocks.*.attn_ln
+		TensorPair attnLn0;
+		// decoder.blocks.*.attn.out
+		TensorPair attnLn1;
+		// decoder.blocks.*.attn.query
+		TensorPair attnQuery;
+		// decoder.blocks.*.attn.key
+		Tensor attnKey;
+		// decoder.blocks.*.attn.value
+		TensorPair attnValue;
+		// decoder.blocks.*.cross_attn_ln
+		TensorPair crossAttnLn0;
+		// decoder.blocks.*.cross_attn.out
+		TensorPair crossAttnLn1;
+		// decoder.blocks.*.cross_attn.query
+		TensorPair crossAttnQuery;
+
+		// decoder.blocks.*.cross_attn.key
+		// Tensor crossAttnKey;
+		// decoder.blocks.*.cross_attn.value
+		// TensorPair crossAttnValue;
+
+		// decoder.blocks.*.mlp_ln
+		TensorPair mlpLn;
+		// decoder.blocks.*.mlp.0
+		TensorPair mlp0;
+		// decoder.blocks.*.mlp.2
+		TensorPair mlp1;
+
+#if TENSOR_GGML_COMPAT
+		// decoder.blocks.*.attn_ln
+		ggml_tensor* attn_ln_0_w;
+		ggml_tensor* attn_ln_0_b;
+
+		// decoder.blocks.*.attn.out
+		ggml_tensor* attn_ln_1_w;
+		ggml_tensor* attn_ln_1_b;
+
+		// decoder.blocks.*.attn.query
+		ggml_tensor* attn_q_w;
+		ggml_tensor* attn_q_b;
+
+		// decoder.blocks.*.attn.key
+		ggml_tensor* attn_k_w;
+
+		// decoder.blocks.*.attn.value
+		ggml_tensor* attn_v_w;
+		ggml_tensor* attn_v_b;
+
+		// decoder.blocks.*.cross_attn_ln
+		ggml_tensor* cross_attn_ln_0_w;
+		ggml_tensor* cross_attn_ln_0_b;
+
+		// decoder.blocks.*.cross_attn.out
+		ggml_tensor* cross_attn_ln_1_w;
+		ggml_tensor* cross_attn_ln_1_b;
+
+		// decoder.blocks.*.cross_attn.query
+		ggml_tensor* cross_attn_q_w;
+		ggml_tensor* cross_attn_q_b;
+
+		// decoder.blocks.*.mlp_ln
+		ggml_tensor* mlp_ln_w;
+		ggml_tensor* mlp_ln_b;
+
+		// decoder.blocks.*.mlp.0
+		ggml_tensor* mlp_0_w;
+		ggml_tensor* mlp_0_b;
+
+		// decoder.blocks.*.mlp.2
+		ggml_tensor* mlp_1_w;
+		ggml_tensor* mlp_1_b;
+#endif
+	};
+
+	struct DecoderTensors
+	{
+		// decoder.positional_embedding
+		Tensor positionalEmbedding;
+
+		// decoder.token_embedding
+		Tensor tokenEmbedding;
+
+		// decoder.ln
+		TensorPair ln;
+		// A vector of layers
+		std::vector<LayerDecoder> layers;
+
+		void setMemoryBuffer( LargeBuffer&& mem ) noexcept
+		{
+			memory = std::move( mem );
+#if TENSOR_GGML_COMPAT
+			makeCompatTensors();
+#endif
+		}
+
+#if TENSOR_GGML_COMPAT
+		void makeCompatTensors();
+
+		// decoder.positional_embedding
+		ggml_tensor* d_pe; // DD
+
+		// decoder.token_embedding
+		ggml_tensor* d_te; // DD
+
+		// decoder.ln
+		ggml_tensor* d_ln_w; // DD
+		ggml_tensor* d_ln_b; // DD
+#endif
+
+	private:
+		// A smart pointer which owns the memory for all the above tensors
+		LargeBuffer memory;
+#if TENSOR_GGML_COMPAT
+		std::vector<ggml_tensor> ggml;
+#endif
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/HybridLoader.cpp b/Whisper/CPU/HybridLoader.cpp
new file mode 100644
index 0000000..e96e3ac
--- /dev/null
+++ b/Whisper/CPU/HybridLoader.cpp
@@ -0,0 +1,140 @@
+#include "stdafx.h"
+#include "HybridLoader.h"
+using namespace CpuCompute;
+using namespace ComLight;
+
+static void populateDecodeTensorsMap( CAtlMap<CStringA, Tensor*>& map, int layersDec, DecoderTensors& dec )
+{
+	dec.layers.resize( layersDec );
+
+	map[ "decoder.positional_embedding" ] = &dec.positionalEmbedding;
+	map[ "decoder.token_embedding.weight" ] = &dec.tokenEmbedding;
+	map[ "decoder.ln.weight" ] = &dec.ln.w;
+	map[ "decoder.ln.bias" ] = &dec.ln.b;
+
+	CStringA tempString;
+	auto add = [ & ]( const char* name, int i, Tensor& t )
+	{
+		tempString.Format( "decoder.blocks.%i.%s", i, name );
+		map[ tempString ] = &t;
+	};
+
+	auto add2 = [ & ]( const char* name, int i, TensorPair& tensors )
+	{
+		tempString.Format( "decoder.blocks.%i.%s.weight", i, name );
+		map[ tempString ] = &tensors.w;
+		tempString.Format( "decoder.blocks.%i.%s.bias", i, name );
+		map[ tempString ] = &tensors.b;
+	};
+
+	for( int i = 0; i < layersDec; i++ )
+	{
+		auto& gpu = dec.layers[ i ];
+		add2( "mlp_ln", i, gpu.mlpLn );
+		add2( "mlp.0", i, gpu.mlp0 );
+		add2( "mlp.2", i, gpu.mlp1 );
+		add2( "attn_ln", i, gpu.attnLn0 );
+		add2( "attn.query", i, gpu.attnQuery );
+		add( "attn.key.weight", i, gpu.attnKey );
+
+		add2( "attn.value", i, gpu.attnValue );
+		add2( "attn.out", i, gpu.attnLn1 );
+
+		add2( "cross_attn_ln", i, gpu.crossAttnLn0 );
+		add2( "cross_attn.query", i, gpu.crossAttnQuery );
+
+		// These 3 tensors are used by the encode() method, to compute cross-attention buffers
+		// Need them in VRAM even for the hybrid model
+		// add( "cross_attn.key.weight", i, gpu.cross_attn_k_w );
+		// add2( "cross_attn.value", i, gpu.cross_attn_v_w, gpu.cross_attn_v_b );
+		add2( "cross_attn.out", i, gpu.crossAttnLn1 );
+	}
+}
+
+HybridLoader::HybridLoader( DecoderTensors& m, int countLayers ) :
+	destination( m )
+{
+	populateDecodeTensorsMap( map, countLayers, destination );
+	pending.reserve( map.GetCount() );
+}
+
+HRESULT HybridLoader::setupTensor( const CStringA& name, int n_dims, int ftype, const std::array<int, 4>& ne, ComLight::iReadStream* stream, int64_t& postponedBytes )
+{
+	auto p = map.Lookup( name );
+	if( nullptr == p )
+		return S_FALSE;
+
+	Tensor& rdi = *p->m_value;
+	PendingTensor& pt = pending.emplace_back();
+
+	__m128i vec = load16( ne.data() );
+	vec = _mm_insert_epi32( vec, 1, 3 );
+	store16( &rdi.ne, vec );
+	rdi.setDenseStrides();
+
+	pt.destPointer = p->m_value;
+	CHECK( stream->getPosition( pt.streamOffset ) );
+	pt.bufferOffset = bufferBytes;
+
+	size_t cbElement;
+	if( ftype == 0 )
+	{
+		rdi.setType( eDataType::FP32 );
+		cbElement = 4;
+	}
+	else
+	{
+		rdi.setType( eDataType::FP16 );
+		cbElement = 2;
+	}
+
+	const size_t totalElts = (size_t)(uint32_t)ne[ 0 ] * (uint32_t)ne[ 1 ] * (uint32_t)ne[ 2 ];
+	if( totalElts * cbElement > UINT_MAX )
+		return DISP_E_OVERFLOW;
+
+	size_t payloadBytes = cbElement * totalElts;
+	pt.payloadBytes = payloadBytes;
+	CHECK( stream->seek( payloadBytes, eSeekOrigin::Current ) );
+	postponedBytes += (int64_t)payloadBytes;
+
+	payloadBytes = ( payloadBytes + 31 ) & ( ~( (size_t)31 ) );
+	bufferBytes += payloadBytes;
+	return S_OK;
+}
+
+HRESULT HybridLoader::completeLoad( ComLight::iReadStream* stream, iLoaderProgressSink& progressSink )
+{
+	if( pending.size() != map.GetCount() )
+	{
+		logError( u8"Not all tensors loaded from model file - expected %zu, got %zu", map.GetCount(), pending.size() );
+		return E_INVALIDARG;
+	}
+
+	LargeBuffer buffer;
+	CHECK( buffer.allocate( bufferBytes ) );
+
+	uint8_t* rdi = buffer.pointer();
+
+	for( const auto& pt : pending )
+	{
+		if( pt.payloadBytes > INT_MAX )
+			return DISP_E_OVERFLOW;
+		CHECK( stream->seek( pt.streamOffset, eSeekOrigin::Begin ) );
+
+		int written = 0;
+		CHECK( stream->read( rdi, (int)pt.payloadBytes, written ) );
+		CHECK( progressSink.gotBytes( (int64_t)pt.payloadBytes ) );
+
+		pt.destPointer->setDataPointer( rdi );
+
+		const size_t cb = ( pt.payloadBytes + 31 ) & ( ~( (size_t)31 ) );
+		rdi += cb;
+	}
+
+	CHECK( buffer.setReadOnly( bufferBytes ) );
+	destination.setMemoryBuffer( std::move( buffer ) );
+
+	constexpr double mulMb = 1.0 / ( 1 << 20 );
+	logDebug( u8"Loaded %zu decoder tensors, %g MB RAM", pending.size(), mulMb * (double)(int64_t)bufferBytes );
+	return S_OK;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/HybridLoader.h b/Whisper/CPU/HybridLoader.h
new file mode 100644
index 0000000..8b12804
--- /dev/null
+++ b/Whisper/CPU/HybridLoader.h
@@ -0,0 +1,37 @@
+#pragma once
+#include "DecoderTensors.h"
+#include <atlstr.h>
+#include <atlcoll.h>
+#include "../../ComLightLib/streams.h"
+
+namespace CpuCompute
+{
+	__interface iLoaderProgressSink
+	{
+		HRESULT gotBytes( int64_t cb );
+	};
+
+	class HybridLoader
+	{
+		DecoderTensors& destination;
+		CAtlMap<CStringA, Tensor*> map;
+		size_t bufferBytes = 0;
+
+		struct alignas( 32 ) PendingTensor
+		{
+			Tensor* destPointer = nullptr;
+			int64_t streamOffset = 0;
+			size_t bufferOffset = 0;
+			size_t payloadBytes = 0;
+		};
+		std::vector<PendingTensor> pending;
+
+	public:
+
+		HybridLoader( DecoderTensors& m, int countLayers );
+
+		HRESULT setupTensor( const CStringA& name, int n_dims, int ftype, const std::array<int, 4>& ne, ComLight::iReadStream* stream, int64_t& postponedBytes );
+
+		HRESULT completeLoad( ComLight::iReadStream* stream, iLoaderProgressSink& progressSink );
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/KvTensors.h b/Whisper/CPU/KvTensors.h
new file mode 100644
index 0000000..a7897d3
--- /dev/null
+++ b/Whisper/CPU/KvTensors.h
@@ -0,0 +1,36 @@
+#pragma once
+#include "Tensor.h"
+#include "LargeBuffer.h"
+#include "../Whisper/sModelParams.h"
+
+namespace CpuCompute
+{
+	class KvTensors
+	{
+		uint16_t* keys = nullptr;
+		uint16_t* values = nullptr;
+		uint32_t size = 0;
+
+		CpuCompute::LargeBuffer memory;
+
+	public:
+		// Create these two large tensors, FP16 precision
+		HRESULT create( const Whisper::sModelParams& mp );
+
+		// A slice of model.memory_cross_k tensor
+		Tensor keysView( uint32_t len, uint32_t off ) const
+		{
+			if( len + off <= size )
+				return Tensor::fromData( keys + off, eDataType::FP16, len );
+			throw E_BOUNDS;
+		}
+
+		// A slice of model.memory_cross_v tensor
+		Tensor valuesView( uint32_t len, uint32_t off ) const
+		{
+			if( len + off <= size )
+				return Tensor::fromData( values + off, eDataType::FP16, len );
+			throw E_BOUNDS;
+		}
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/KvTensorsCpu.cpp b/Whisper/CPU/KvTensorsCpu.cpp
new file mode 100644
index 0000000..88a70ff
--- /dev/null
+++ b/Whisper/CPU/KvTensorsCpu.cpp
@@ -0,0 +1,19 @@
+#include "stdafx.h"
+#include "KvTensors.h"
+using namespace CpuCompute;
+
+// Create these two large tensors, FP16 precision
+HRESULT KvTensors::create( const Whisper::sModelParams& mp )
+{
+	const uint32_t n_mem = mp.n_text_layer * mp.n_text_ctx;
+	const uint32_t n_elements = mp.n_text_state * n_mem;
+
+	const size_t cb = sizeof( uint16_t ) * (size_t)n_elements * 2;
+	CHECK( memory.allocate( cb ) );
+
+	uint16_t* pointer = (uint16_t*)memory.pointer();
+	keys = pointer;
+	values = pointer + n_elements;
+	size = n_elements;
+	return S_OK;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/LargeBuffer.cpp b/Whisper/CPU/LargeBuffer.cpp
new file mode 100644
index 0000000..e124686
--- /dev/null
+++ b/Whisper/CPU/LargeBuffer.cpp
@@ -0,0 +1,34 @@
+#include "stdafx.h"
+#include "LargeBuffer.h"
+using namespace CpuCompute;
+
+void LargeBuffer::deallocate()
+{
+	if( nullptr == pv )
+		return;
+	VirtualFree( pv, 0, MEM_RELEASE );
+	pv = nullptr;
+}
+
+HRESULT LargeBuffer::allocate( size_t cb )
+{
+	deallocate();
+	
+	pv = VirtualAlloc( nullptr, cb, MEM_COMMIT | MEM_RESERVE, PAGE_READWRITE );
+	if( nullptr != pv )
+		return S_OK;
+	return HRESULT_FROM_WIN32( GetLastError() );
+}
+
+HRESULT LargeBuffer::setReadOnly( size_t cb )
+{
+	if( nullptr != pv )
+	{
+		DWORD op = 0;
+		if( VirtualProtect( pv, cb, PAGE_READONLY, &op ) )
+			return S_OK;
+		return HRESULT_FROM_WIN32( GetLastError() );
+	}
+	else
+		return OLE_E_BLANK;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/LargeBuffer.h b/Whisper/CPU/LargeBuffer.h
new file mode 100644
index 0000000..d9a9a22
--- /dev/null
+++ b/Whisper/CPU/LargeBuffer.h
@@ -0,0 +1,44 @@
+#pragma once
+
+namespace CpuCompute
+{
+	// A large memory buffer allocated with VirtualAlloc kernel API, bypassing the heap.
+	class LargeBuffer
+	{
+		void* pv = nullptr;
+	public:
+		LargeBuffer() = default;
+		LargeBuffer( const LargeBuffer& ) = delete;
+		LargeBuffer( LargeBuffer&& that ) noexcept
+		{
+			pv = that.pv;
+			that.pv = nullptr;
+		}
+		~LargeBuffer()
+		{
+			deallocate();
+		}
+		void operator=( LargeBuffer&& that ) noexcept
+		{
+			std::swap( pv, that.pv );
+		}
+		void operator=( const LargeBuffer& that ) = delete;
+
+		// Allocate buffer with specified count of bytes, and read+write memory protection
+		// The OS kernel guarantees zero-initialization of that memory.
+		HRESULT allocate( size_t cb );
+
+		// Change memory protection of the buffer to read only
+		HRESULT setReadOnly( size_t cb );
+
+		// Unless the pointer is nullptr, deallocate the buffer
+		void deallocate();
+
+		// Pointer to the start of the buffer, aligned by memory page = 4 kilobytes
+		uint8_t* pointer() const
+		{
+			assert( nullptr != pv );
+			return (uint8_t*)pv;
+		}
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/MlContext.h b/Whisper/CPU/MlContext.h
new file mode 100644
index 0000000..062a7b8
--- /dev/null
+++ b/Whisper/CPU/MlContext.h
@@ -0,0 +1,71 @@
+#pragma once
+#include "Tensor.h"
+#include "ParallelForRunner.h"
+
+namespace CpuCompute
+{
+	class MlContext
+	{
+		ParallelForRunner pfor;
+		iMemoryAllocator* allocator = nullptr;
+
+	public:
+		MlContext( int threads );
+		MlContext( const MlContext& ) = delete;
+		~MlContext() = default;
+
+		HRESULT setThreadsCount( int threads )
+		{
+			return pfor.setThreadsCount( threads );
+		}
+
+		iMemoryAllocator* setAllocator( iMemoryAllocator* alloc )
+		{
+			iMemoryAllocator* const ret = allocator;
+			allocator = alloc;
+			return ret;
+		}
+
+		Tensor createTensor( eDataType type, const std::array<uint32_t, 4>& size );
+		Tensor createTensor( eDataType type, std::initializer_list<uint32_t> size );
+
+		Tensor addRows( const Tensor& d_te, const Tensor& d_pe, const int* tokens, const int n_tokens, const int n_past );
+
+		Tensor norm( const Tensor& arg );
+
+		// cur = add( mul( repeat( w, cur ), cur ), repeat( b, cur ) );
+		void fmaRepeat( Tensor& cur, const Tensor& w, const Tensor& b );
+
+		inline void fmaRepeat( Tensor& cur, const TensorPair wb )
+		{
+			fmaRepeat( cur, wb.w, wb.b );
+		}
+
+		// Multiply two matrices
+		Tensor mulMat( const Tensor& a, const Tensor& b );
+
+		// cur = add( repeat( b, cur ), cur ); cur = scale(cur, scaling)
+		void addRepeatScale( Tensor& cur, const Tensor& b, float scaling );
+
+		void addRepeat( Tensor& cur, const Tensor& b );
+
+		Tensor add( const Tensor& a, const Tensor& b );
+		void addInPlace( Tensor& a, const Tensor& b );
+		void addRepeatGelu( Tensor& cur, const Tensor& b );
+
+		// cur = scale(cur, scaling)
+		void scale( Tensor& cur, float scaling );
+
+		void diagMaskInf( Tensor& cur, uint32_t n_past );
+
+		void softMax( Tensor& cur, float inputScale = 1.0f );
+
+		Tensor copy( const Tensor& a, eDataType type, std::initializer_list<uint32_t> size );
+
+		HRESULT copyImpl( Tensor& result, const Tensor& source );
+
+		Tensor permute( const Tensor& a, uint8_t axis0, uint8_t axis1, uint8_t axis2, uint8_t axis3 );
+
+		void copyInPlace( Tensor& dest, const Tensor& a, eDataType type, std::initializer_list<uint32_t> size );
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/MlContextCpu.cpp b/Whisper/CPU/MlContextCpu.cpp
new file mode 100644
index 0000000..e34a823
--- /dev/null
+++ b/Whisper/CPU/MlContextCpu.cpp
@@ -0,0 +1,597 @@
+#include "stdafx.h"
+#include "MlContext.h"
+#include "simdUtils.h"
+#include "mulMat.h"
+using namespace CpuCompute;
+
+MlContext::MlContext( int threads ) : pfor( threads )
+{
+}
+
+Tensor MlContext::createTensor( eDataType type, const std::array<uint32_t, 4>& size )
+{
+	Tensor res;
+	check( res.create( type, size, allocator ) );
+	return res;
+}
+
+Tensor MlContext::createTensor( eDataType type, std::initializer_list<uint32_t> size )
+{
+	Tensor res;
+	check( res.create( type, size, allocator ) );
+	return res;
+}
+
+namespace
+{
+	inline const uint16_t* getRow16( const Tensor& t, size_t index )
+	{
+		const uint16_t* rsi = t.fp16();
+		rsi += index * t.nb[ 1 ];
+		return rsi;
+	}
+	inline const float* getRow32( const Tensor& t, size_t index )
+	{
+		const float* rsi = t.fp32();
+		rsi += index * t.nb[ 1 ];
+		return rsi;
+	}
+}
+
+Tensor MlContext::addRows( const Tensor& d_te, const Tensor& d_pe, const int* tokens, const int n_tokens, const int n_past )
+{
+	if( d_te.type() != eDataType::FP16 || d_pe.type() != eDataType::FP32 )
+		throw E_INVALIDARG;
+	if( d_te.ne[ 0 ] != d_pe.ne[ 0 ] )
+		throw E_INVALIDARG;
+	if( n_tokens <= 0 )
+		throw E_BOUNDS;
+
+	Tensor res = createTensor( eDataType::FP32, { d_te.ne[ 0 ], (uint32_t)n_tokens } );
+
+	const size_t inner = (size_t)d_te.ne[ 0 ];
+	const size_t outer = (size_t)n_tokens;
+	float* rdi = res.fp32();
+	for( size_t i = 0; i < outer; i++, rdi += inner, tokens++ )
+	{
+		const uint16_t* const source1 = getRow16( d_te, *(const uint32_t*)tokens );
+		const float* const source2 = getRow32( d_pe, i + (size_t)n_past );
+		addF16to32( rdi, source1, source2, inner );
+	}
+	return res;
+}
+
+namespace
+{
+	class DispatchHelper3
+	{
+		std::array<uint32_t, 3> ne;
+
+	public:
+		DispatchHelper3() = default;
+		DispatchHelper3( uint32_t x, uint32_t y, uint32_t z )
+		{
+			assert( x > 0 && y > 0 && z > 0 );
+			ne[ 0 ] = x;
+			ne[ 1 ] = y;
+			ne[ 2 ] = z;
+		}
+		size_t groupsCount() const
+		{
+			size_t res = ne[ 0 ];
+			res *= ne[ 1 ];
+			res *= ne[ 2 ];
+			return res;
+		}
+		std::array<uint32_t, 3> unpack( size_t idx ) const
+		{
+			assert( idx < groupsCount() );
+			std::array<uint32_t, 3> res;
+			res[ 0 ] = (uint32_t)( idx % ne[ 0 ] );
+			idx = idx / ne[ 0 ];
+			res[ 1 ] = (uint32_t)( idx % ne[ 1 ] );
+			res[ 2 ] = (uint32_t)( idx / ne[ 1 ] );
+			return res;
+		}
+		void next( std::array<uint32_t, 3>& i ) const
+		{
+			i[ 0 ]++;
+			if( i[ 0 ] < ne[ 0 ] )
+				return;
+			i[ 0 ] = 0;
+			i[ 1 ]++;
+			if( i[ 1 ] < ne[ 1 ] )
+				return;
+			i[ 1 ] = 0;
+			i[ 2 ]++;
+		}
+	};
+
+	inline const float* sourceRow( const float* rsi, const std::array<uint32_t, 3>& idx, size_t nb0, size_t nb1, size_t nb2 )
+	{
+		const size_t r0 = idx[ 0 ] * nb0;
+		const size_t r1 = idx[ 1 ] * nb1;
+		const size_t r2 = idx[ 2 ] * nb2;
+		rsi = rsi + r0 + r1 + r2;
+		return rsi;
+	}
+
+	struct NormContext : public iComputeRange
+	{
+		const float* source;
+		float* result;
+		size_t inner;
+		DispatchHelper3 threads;
+		std::array<uint32_t, 3> nbInput;
+
+		HRESULT __stdcall compute( size_t i, size_t end ) const override final
+		{
+			ALIGNED_SPAN( temp, inner );
+
+			std::array<uint32_t, 3> idx = threads.unpack( i );
+			float* rdi = result + i * inner;
+			for( ; i < end; i++, rdi += inner, threads.next( idx ) )
+			{
+				const float* rsi = sourceRow( source, idx, nbInput[ 0 ], nbInput[ 1 ], nbInput[ 2 ] );
+				norm( rdi, temp, rsi, inner );
+			}
+			return S_OK;
+		}
+	};
+}
+
+Tensor MlContext::norm( const Tensor& arg )
+{
+	if( arg.type() != eDataType::FP32 || arg.nb[ 0 ] != 1 )
+		throw E_INVALIDARG;
+	Tensor res = createTensor( eDataType::FP32, arg.ne );
+
+	NormContext context;
+	context.source = arg.fp32();
+	context.result = res.fp32();
+	context.inner = arg.ne[ 0 ];
+	context.threads = DispatchHelper3( arg.ne[ 1 ], arg.ne[ 2 ], arg.ne[ 3 ] );
+	context.nbInput = { arg.nb[ 1 ], arg.nb[ 2 ], arg.nb[ 3 ] };
+
+	check( pfor.parallelFor( context, context.threads.groupsCount() ) );
+	return res;
+}
+
+void MlContext::fmaRepeat( Tensor& cur, const Tensor& w, const Tensor& b )
+{
+	if( !( cur.isContinuous() && w.isContinuous() && b.isContinuous() ) )
+		throw E_INVALIDARG;
+
+	if( !( cur.type() == eDataType::FP32 && w.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	if( !isSameShape( w, b ) )
+		throw E_INVALIDARG;
+
+	DispatchHelper3 helper{ cur.ne[ 1 ], cur.ne[ 2 ], cur.ne[ 3 ] };
+	std::array<uint32_t, 3> idx = { 0, 0, 0 };
+	const size_t countRows = helper.groupsCount();
+
+	const size_t innerRes = cur.ne[ 0 ];
+	const size_t innerPattern = w.ne[ 0 ];
+
+	float* rdi = cur.fp32();
+	for( size_t i = 0; i < countRows; i++, helper.next( idx ), rdi += innerRes )
+	{
+		std::array<uint32_t, 3> idxPattern;
+		idxPattern[ 0 ] = idx[ 0 ] % w.ne[ 1 ];
+		idxPattern[ 1 ] = idx[ 1 ] % w.ne[ 2 ];
+		idxPattern[ 2 ] = idx[ 2 ] % w.ne[ 3 ];
+
+		const float* s1 = sourceRow( w.fp32(), idxPattern, w.nb[ 1 ], w.nb[ 2 ], w.nb[ 3 ] );
+		const float* s2 = sourceRow( b.fp32(), idxPattern, b.nb[ 1 ], b.nb[ 2 ], b.nb[ 3 ] );
+		fmaRepeatRow( rdi, innerRes, s1, s2, innerPattern );
+	}
+}
+
+Tensor MlContext::mulMat( const Tensor& a, const Tensor& b )
+{
+	if( !DirectCompute::canMulMat( a, b ) )
+		throw E_INVALIDARG;
+
+	std::array<uint32_t, 4> ne{ a.ne[ 1 ], b.ne[ 1 ], a.ne[ 2 ], b.ne[ 3 ] };
+	Tensor result = createTensor( eDataType::FP32, ne );
+
+	check( CpuCompute::mulMat( result, a, b, pfor ) );
+	return result;
+}
+
+// cur = add( repeat( b, cur ), cur ); cur = scale(cur, scaling)
+void MlContext::addRepeatScale( Tensor& cur, const Tensor& b, float scaling )
+{
+	if( !( cur.isContinuous() && b.isContinuous() ) )
+		throw E_INVALIDARG;
+	if( !( cur.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	DispatchHelper3 helper{ cur.ne[ 1 ], cur.ne[ 2 ], cur.ne[ 3 ] };
+	std::array<uint32_t, 3> idx = { 0, 0, 0 };
+	const size_t countRows = helper.groupsCount();
+
+	const size_t innerRes = (uint32_t)cur.ne[ 0 ];
+	const size_t innerPattern = (uint32_t)b.ne[ 0 ];
+
+	float* rdi = cur.fp32();
+	const __m256 scale = _mm256_set1_ps( scaling );
+	for( size_t i = 0; i < countRows; i++, helper.next( idx ), rdi += innerRes )
+	{
+		std::array<uint32_t, 3> idxPattern;
+		idxPattern[ 0 ] = idx[ 0 ] % (uint32_t)b.ne[ 1 ];
+		idxPattern[ 1 ] = idx[ 1 ] % (uint32_t)b.ne[ 2 ];
+		idxPattern[ 2 ] = idx[ 2 ] % (uint32_t)b.ne[ 3 ];
+
+		const float* source = sourceRow( b.fp32(), idxPattern, b.nb[ 1 ], b.nb[ 2 ], b.nb[ 3 ] );
+		addRepeatScaleRow( rdi, innerRes, source, innerPattern, scale );
+	}
+}
+
+void MlContext::addRepeat( Tensor& cur, const Tensor& b )
+{
+	if( !( cur.isContinuous() && b.isContinuous() ) )
+		throw E_INVALIDARG;
+	if( !( cur.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	DispatchHelper3 helper{ cur.ne[ 1 ], cur.ne[ 2 ], cur.ne[ 3 ] };
+	std::array<uint32_t, 3> idx = { 0, 0, 0 };
+	const size_t countRows = helper.groupsCount();
+
+	const size_t innerRes = (uint32_t)cur.ne[ 0 ];
+	const size_t innerPattern = (uint32_t)b.ne[ 0 ];
+
+	float* rdi = cur.fp32();
+	for( size_t i = 0; i < countRows; i++, helper.next( idx ), rdi += innerRes )
+	{
+		std::array<uint32_t, 3> idxPattern;
+		idxPattern[ 0 ] = idx[ 0 ] % (uint32_t)b.ne[ 1 ];
+		idxPattern[ 1 ] = idx[ 1 ] % (uint32_t)b.ne[ 2 ];
+		idxPattern[ 2 ] = idx[ 2 ] % (uint32_t)b.ne[ 3 ];
+
+		const float* source = sourceRow( b.fp32(), idxPattern, b.nb[ 1 ], b.nb[ 2 ], b.nb[ 3 ] );
+		addRepeatRow( rdi, innerRes, source, innerPattern );
+	}
+}
+
+// cur = scale(cur, scaling)
+void MlContext::scale( Tensor& cur, float scaling )
+{
+	if( !( cur.isContinuous() && cur.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	const size_t len = cur.countElements();
+	const __m256 scale = _mm256_set1_ps( scaling );
+	scaleRow( cur.fp32(), len, scale );
+}
+
+void MlContext::diagMaskInf( Tensor& cur, uint32_t n_past )
+{
+	if( !( cur.isContinuous() && cur.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	const size_t n = cur.countRows();
+	const size_t nc = cur.ne[ 0 ];
+	const size_t nr = cur.ne[ 1 ];
+	const size_t nz = n / nr;
+
+	for( size_t k = 0; k < nz; k++ )
+	{
+		for( size_t j = 0; j < nr; j++ )
+		{
+			float* const rdi = cur.fp32() + k * cur.nb[ 2 ] + j * cur.nb[ 1 ];
+			// +1 because the original code checked for `if( i > n_past + j )`
+			// That's why the first index to write is ( n_past + j + 1 )
+			const size_t start = n_past + j + 1;
+			const ptrdiff_t len = (ptrdiff_t)nc - (ptrdiff_t)start;
+			if( len <= 0 )
+				continue;
+
+			// Generates a store string instruction (rep stosd).
+			// The magic number is negative infinity in FP32: https://www.h-schmidt.net/FloatConverter/IEEE754.html
+			__stosd( (DWORD*)( rdi + start ), 0xff800000u, (size_t)len );
+		}
+	}
+}
+
+void MlContext::softMax( Tensor& cur, float inputScale )
+{
+	if( !( cur.isContinuous() && cur.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	struct SoftMaxContext : public iComputeRange
+	{
+		float* data;
+		float inputScale;
+		size_t length, stride;
+
+		HRESULT __stdcall compute( size_t i, size_t end ) const override final
+		{
+			float* rdi = data + stride * i;
+			for( ; i < end; i++, rdi += stride )
+				::softMax( rdi, length, inputScale );
+			return S_OK;
+		}
+	};
+
+	SoftMaxContext context;
+	context.data = cur.fp32();
+	context.inputScale = inputScale;
+	context.length = cur.ne[ 0 ];
+	context.stride = cur.nb[ 1 ];
+
+	const size_t n = cur.countRows();
+	pfor.parallelFor( context, n );
+}
+
+namespace
+{
+	template<class R, class S>
+	__forceinline void copyElement( R* rdi, const S* rsi )
+	{
+		static_assert( std::is_same<R, S>() );
+		*rdi = *rsi;
+	}
+	template<>
+	__forceinline void copyElement<float, uint16_t>( float* rdi, const uint16_t* rsi )
+	{
+		__m128i iv = _mm_cvtsi32_si128( *rsi );
+		__m128 fv = _mm_cvtph_ps( iv );
+		_mm_store_ss( rdi, fv );
+	}
+	template<>
+	__forceinline void copyElement<uint16_t, float>( uint16_t* rdi, const float* rsi )
+	{
+		__m128 fv = _mm_load_ss( rsi );
+		__m128i iv = _mm_cvtps_ph( fv, 0 );
+		*rdi = (uint16_t)(uint32_t)_mm_cvtsi128_si32( iv );
+	}
+
+	template<class R, class S>
+	__forceinline void copyRow( R* rdi, const S* rsi, size_t length )
+	{
+		static_assert( std::is_same<R, S>() );
+		memcpy( rdi, rsi, length * sizeof( R ) );
+	}
+	template<>
+	__forceinline void copyRow<uint16_t, float>( uint16_t* rdi, const float* rsi, size_t length )
+	{
+		floatsDowncast( rdi, rsi, length );
+	}
+	template<>
+	__forceinline void copyRow<float, uint16_t>( float* rdi, const uint16_t* rsi, size_t length )
+	{
+		floatsUpcast( rdi, rsi, length );
+	}
+
+	template<class R, class S>
+	static void __declspec( noinline ) copyImpl( R* rdi, const S* rsi, const TensorShape& shape )
+	{
+		const bool continuousRows = shape.nb[ 0 ] == 1;
+
+		for( size_t i03 = 0; i03 < shape.ne[ 3 ]; i03++, rsi += shape.nb[ 3 ] )
+		{
+			const S* source2 = rsi;
+			for( size_t i02 = 0; i02 < shape.ne[ 2 ]; i02++, source2 += shape.nb[ 2 ] )
+			{
+				const S* source1 = source2;
+				for( size_t i01 = 0; i01 < shape.ne[ 1 ]; i01++, source1 += shape.nb[ 1 ] )
+				{
+					// Performance optimization here: when the rows are dense, we can copy them much faster with memcpy()
+					// Or at least with AVX, when we need to convert between numeric types
+					if( continuousRows )
+					{
+						// This branch is very predictable, same outcome for all loop iterations
+						copyRow( rdi, source1, shape.ne[ 0 ] );
+						rdi += shape.ne[ 0 ];
+					}
+					else
+					{
+						const S* source0 = source1;
+						for( size_t i00 = 0; i00 < shape.ne[ 0 ]; i00++, source0 += shape.nb[ 0 ] )
+						{
+							copyElement( rdi, source0 );
+							rdi++;
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+HRESULT MlContext::copyImpl( Tensor& result, const Tensor& source )
+{
+	if( !( result.isContinuous() && ( result.countElements() == source.countElements() ) ) )
+		return E_INVALIDARG;
+
+	const eDataType typeResult = result.type();
+	const eDataType typeSource = source.type();
+	if( source.isContinuous() )
+	{
+		const size_t elts = result.countElements();
+		if( typeResult == typeSource )
+		{
+			const size_t bytes = elts * elementSize( typeResult );
+			memcpy( result.data(), source.data(), bytes );
+			return S_OK;
+		}
+		if( typeSource == eDataType::FP16 && typeResult == eDataType::FP32 )
+		{
+			floatsUpcast( result.fp32(), source.fp16(), elts );
+			return S_OK;
+		}
+		if( typeSource == eDataType::FP32 && typeResult == eDataType::FP16 )
+		{
+			floatsDowncast( result.fp16(), source.fp32(), elts );
+			return S_OK;
+		}
+		return E_UNEXPECTED;
+	}
+	else
+	{
+		if( typeSource == eDataType::FP16 && typeResult == eDataType::FP16 )
+		{
+			::copyImpl( result.fp16(), source.fp16(), source );
+			return S_OK;
+		}
+		if( typeSource == eDataType::FP32 && typeResult == eDataType::FP32 )
+		{
+			::copyImpl( result.fp32(), source.fp32(), source );
+			return S_OK;
+		}
+		if( typeSource == eDataType::FP16 && typeResult == eDataType::FP32 )
+		{
+			::copyImpl( result.fp32(), source.fp16(), source );
+			return S_OK;
+		}
+		if( typeSource == eDataType::FP32 && typeResult == eDataType::FP16 )
+		{
+			::copyImpl( result.fp16(), source.fp32(), source );
+			return S_OK;
+		}
+		return E_UNEXPECTED;
+	}
+}
+
+Tensor MlContext::copy( const Tensor& a, eDataType type, std::initializer_list<uint32_t> size )
+{
+	const size_t dims = size.size();
+	if( 0 == dims || dims > 4 )
+		throw E_BOUNDS;
+
+	size_t nRequested = 1;
+	for( size_t i = 0; i < dims; i++ )
+	{
+		uint32_t n = size.begin()[ i ];
+		nRequested *= n;
+	}
+	if( nRequested != a.countElements() )
+		throw E_INVALIDARG;
+
+	if( a.type() == type && a.isContinuous() )
+	{
+		// Same type, and it's dense - no need to move data, equal to reshape
+		Tensor res{ a };
+		for( size_t i = 0; i < dims; i++ )
+			res.ne[ i ] = size.begin()[ i ];;
+		for( size_t i = dims; i < 4; i++ )
+			res.ne[ i ] = 1;
+		res.setDenseStrides();
+		return res;
+	}
+	else
+	{
+		// Need to convert types, and/or transpose the tensor. Make another tensor for the output
+		Tensor res = createTensor( type, size );
+		check( copyImpl( res, a ) );
+		return res;
+	}
+}
+
+Tensor MlContext::permute( const Tensor& a, uint8_t axis0, uint8_t axis1, uint8_t axis2, uint8_t axis3 )
+{
+	assert( axis0 < 4 );
+	assert( axis1 < 4 );
+	assert( axis2 < 4 );
+	assert( axis3 < 4 );
+
+	assert( axis0 != axis1 );
+	assert( axis0 != axis2 );
+	assert( axis0 != axis3 );
+	assert( axis1 != axis2 );
+	assert( axis1 != axis3 );
+	assert( axis2 != axis3 );
+
+	Tensor res = a;
+	res.ne[ axis0 ] = a.ne[ 0 ];
+	res.ne[ axis1 ] = a.ne[ 1 ];
+	res.ne[ axis2 ] = a.ne[ 2 ];
+	res.ne[ axis3 ] = a.ne[ 3 ];
+
+	res.nb[ axis0 ] = a.nb[ 0 ];
+	res.nb[ axis1 ] = a.nb[ 1 ];
+	res.nb[ axis2 ] = a.nb[ 2 ];
+	res.nb[ axis3 ] = a.nb[ 3 ];
+
+	return res;
+}
+
+void MlContext::copyInPlace( Tensor& dest, const Tensor& a, eDataType type, std::initializer_list<uint32_t> size )
+{
+	assert( type == dest.type() );
+
+	const size_t dims = size.size();
+	if( 0 == dims || dims > 4 )
+		throw E_BOUNDS;
+
+	size_t nRequested = 1;
+	for( size_t i = 0; i < dims; i++ )
+	{
+		uint32_t n = size.begin()[ i ];
+		nRequested *= n;
+	}
+	if( nRequested != a.countElements() || nRequested != dest.countElements() )
+		throw E_INVALIDARG;
+
+	// Reshape the destination
+	for( size_t i = 0; i < dims; i++ )
+		dest.ne[ i ] = size.begin()[ i ];
+	for( size_t i = dims; i < 4; i++ )
+		dest.ne[ i ] = 1;
+	dest.setDenseStrides();
+
+	// Copy the data
+	check( copyImpl( dest, a ) );
+}
+
+void MlContext::addInPlace( Tensor& a, const Tensor& b )
+{
+	if( !( a.isContinuous() && b.isContinuous() && a.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_NOTIMPL;
+
+	const size_t length = a.countElements();
+	addRowInPlace( a.fp32(), b.fp32(), length );
+}
+
+Tensor MlContext::add( const Tensor& a, const Tensor& b )
+{
+	if( !( a.isContinuous() && b.isContinuous() && a.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_NOTIMPL;
+
+	Tensor res = createTensor( eDataType::FP32, a.ne );
+	const size_t length = a.countElements();
+	addRow( res.fp32(), a.fp32(), b.fp32(), length );
+	return res;
+}
+
+void MlContext::addRepeatGelu( Tensor& cur, const Tensor& b )
+{
+	if( !( cur.isContinuous() && b.isContinuous() ) )
+		throw E_INVALIDARG;
+	if( !( cur.type() == eDataType::FP32 && b.type() == eDataType::FP32 ) )
+		throw E_INVALIDARG;
+
+	DispatchHelper3 helper{ cur.ne[ 1 ], cur.ne[ 2 ], cur.ne[ 3 ] };
+	std::array<uint32_t, 3> idx = { 0, 0, 0 };
+	const size_t countRows = helper.groupsCount();
+
+	const size_t innerRes = (uint32_t)cur.ne[ 0 ];
+	const size_t innerPattern = (uint32_t)b.ne[ 0 ];
+	float* rdi = cur.fp32();
+	auto& lookupTables = getLookupTables();
+	for( size_t i = 0; i < countRows; i++, helper.next( idx ), rdi += innerRes )
+	{
+		std::array<uint32_t, 3> idxPattern;
+		idxPattern[ 0 ] = idx[ 0 ] % (uint32_t)b.ne[ 1 ];
+		idxPattern[ 1 ] = idx[ 1 ] % (uint32_t)b.ne[ 2 ];
+		idxPattern[ 2 ] = idx[ 2 ] % (uint32_t)b.ne[ 3 ];
+
+		const float* source = sourceRow( b.fp32(), idxPattern, b.nb[ 1 ], b.nb[ 2 ], b.nb[ 3 ] );
+		addRepeatGeluRow( rdi, innerRes, source, innerPattern, lookupTables );
+	}
+	return;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/ParallelForRunner.cpp b/Whisper/CPU/ParallelForRunner.cpp
new file mode 100644
index 0000000..7151a23
--- /dev/null
+++ b/Whisper/CPU/ParallelForRunner.cpp
@@ -0,0 +1,149 @@
+#include "stdafx.h"
+#include "ParallelForRunner.h"
+using namespace CpuCompute;
+
+ParallelForRunner::ParallelForRunner( int threads ) :
+	maxThreads( threads )
+{
+	if( maxThreads <= 1 )
+	{
+		threadBuffers.resize( 1 );
+		return;
+	}
+
+	work = CreateThreadpoolWork( &workCallbackStatic, this, nullptr );
+	if( nullptr == work )
+		throw getLastHr();
+	threadBuffers.resize( maxThreads );
+}
+
+HRESULT ParallelForRunner::setThreadsCount( int threads )
+{
+	maxThreads = threads;
+	if( threads <= 1 )
+	{
+		threadBuffers.resize( 1 );
+		return S_OK;
+	}
+
+	threadBuffers.resize( maxThreads );
+	if( nullptr == work )
+	{
+		work = CreateThreadpoolWork( &workCallbackStatic, this, nullptr );
+		if( nullptr == work )
+			return getLastHr();
+	}
+	return S_OK;
+}
+
+ParallelForRunner::~ParallelForRunner()
+{
+	if( nullptr != work )
+	{
+		if( S_FALSE == status )
+			WaitForThreadpoolWorkCallbacks( work, FALSE );
+		CloseThreadpoolWork( work );
+	}
+}
+
+namespace
+{
+	thread_local uint32_t currentThreadIndex = UINT_MAX;
+}
+
+void ParallelForRunner::runBatch( size_t ith ) noexcept
+{
+	currentThreadIndex = (uint32_t)ith;
+	const size_t begin = ( ith * countItems ) / countThreads;
+	const size_t end = ( ( ith + 1 ) * countItems ) / countThreads;
+
+	HRESULT hr = E_UNEXPECTED;
+	try
+	{
+		hr = computeRange->compute( begin, end );
+	}
+	catch( HRESULT code )
+	{
+		hr = code;
+	}
+	catch( const std::bad_alloc& )
+	{
+		hr = E_OUTOFMEMORY;
+	}
+	catch( const std::exception& )
+	{
+		hr = E_FAIL;
+	}
+	currentThreadIndex = UINT_MAX;
+	if( SUCCEEDED( hr ) )
+		return;
+	InterlockedCompareExchange( &status, hr, S_FALSE );
+}
+
+void* ParallelForRunner::threadLocalBuffer( size_t cb )
+{
+	const uint32_t idx = currentThreadIndex;
+	if( idx < threadBuffers.size() )
+	{
+		ThreadBuffer& tb = threadBuffers[ idx ];
+		if( tb.cb >= cb )
+		{
+			// We already have large enough buffer for the current thread
+			return tb.memory.pointer();
+		}
+		tb.memory.deallocate();
+		check( tb.memory.allocate( cb ) );
+		tb.cb = cb;
+		return tb.memory.pointer();
+	}
+	if( idx != UINT_MAX )
+		throw E_BOUNDS;
+	else
+	{
+		logError( u8"threadLocalBuffer() method only works from inside a pool callback" );
+		throw E_UNEXPECTED;
+	}
+}
+
+void __stdcall ParallelForRunner::workCallbackStatic( PTP_CALLBACK_INSTANCE Instance, void* pv, PTP_WORK Work ) noexcept
+{
+	ParallelForRunner& context = *(ParallelForRunner*)pv;
+	const size_t ith = (uint32_t)( InterlockedIncrement( &context.threadIndex ) );
+	context.runBatch( ith );
+}
+
+HRESULT ParallelForRunner::parallelFor( iComputeRange& compute, size_t length, size_t minBatch )
+{
+	if( maxThreads <= 1 )
+	{
+		currentThreadIndex = 0;
+		const HRESULT hr1 = compute.compute( 0, length );
+		currentThreadIndex = UINT_MAX;
+		return hr1;
+	}
+	assert( minBatch > 0 );
+
+	size_t nth = length / minBatch;
+	nth = std::min( nth, (size_t)(uint32_t)maxThreads );
+
+	computeRange = &compute;
+	countItems = length;
+	countThreads = nth;
+	threadIndex = 0;
+	status = S_FALSE;
+
+	for( size_t i = 1; i < nth; i++ )
+		SubmitThreadpoolWork( work );
+	runBatch( 0 );
+
+	if( nth > 1 )
+		WaitForThreadpoolWorkCallbacks( work, FALSE );
+
+	computeRange = nullptr;
+	const HRESULT hr = status;
+	status = S_OK;
+	if( SUCCEEDED( hr ) )
+		return S_OK;
+
+	return hr;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/ParallelForRunner.h b/Whisper/CPU/ParallelForRunner.h
new file mode 100644
index 0000000..baef647
--- /dev/null
+++ b/Whisper/CPU/ParallelForRunner.h
@@ -0,0 +1,52 @@
+#pragma once
+#include "LargeBuffer.h"
+
+namespace CpuCompute
+{
+	// Callback interface for the parallel `for`
+	__interface iComputeRange
+	{
+		// The implementation calls this method on multiple thread pool threads in parallel, and aggregates status codes.
+		HRESULT __stdcall compute( size_t begin, size_t end ) const;
+	};
+
+	// Similar to ThreadPoolWork in parallelFor.h, optimized to be used as a direct replacement of OpenMP pool.
+	class alignas( 64 ) ParallelForRunner
+	{
+	public:
+		ParallelForRunner( int threads );
+		~ParallelForRunner();
+
+		HRESULT setThreadsCount( int threads );
+
+		HRESULT parallelFor( iComputeRange& compute, size_t length, size_t minBatch = 1 );
+
+		// Allocate a temporary buffer for the calling thread.
+		// The pointer is guaranteed to be aligned by page size = 4kb
+		void* threadLocalBuffer( size_t cb );
+
+	private:
+
+		int maxThreads;
+		PTP_WORK work = nullptr;
+		iComputeRange* computeRange = nullptr;
+		size_t countItems = 0;
+		size_t countThreads = 0;
+
+		// Aligning by cache lines.
+		// Avoiding cache line sharing between CPU cores improves performance, despite wasting a few bytes of memory.
+		struct alignas( 64 ) ThreadBuffer
+		{
+			LargeBuffer memory;
+			size_t cb = 0;
+		};
+		std::vector<ThreadBuffer> threadBuffers;
+
+		alignas( 64 ) volatile long threadIndex = 0;
+		volatile HRESULT status = S_OK;
+
+		void runBatch( size_t ith ) noexcept;
+
+		static void __stdcall workCallbackStatic( PTP_CALLBACK_INSTANCE Instance, void* pv, PTP_WORK Work ) noexcept;
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/Readme.txt b/Whisper/CPU/Readme.txt
new file mode 100644
index 0000000..702813d
--- /dev/null
+++ b/Whisper/CPU/Readme.txt
@@ -0,0 +1 @@
+The code in this folder is dropped by the linker’s dead code elimination optimization pass, unless you change BUILD_HYBRID_VERSION macro in stdafx.h
\ No newline at end of file
diff --git a/Whisper/CPU/Tensor.h b/Whisper/CPU/Tensor.h
new file mode 100644
index 0000000..2fcde63
--- /dev/null
+++ b/Whisper/CPU/Tensor.h
@@ -0,0 +1,139 @@
+#pragma once
+#include "../D3D/enums.h"
+#include "../ML/TensorShape.h"
+// 1 = new tensors can be allocated with `nullptr` iMemoryAllocator, by allocating memory internally and counting these references
+// 0 = memory allocator is mandatory, create() methods will fail with E_POINTER if the allocator is `nullptr`
+#define TENSOR_INTERNAL_ALLOC 0
+
+// 1 = expose compatibility API for GGML interop
+#define TENSOR_GGML_COMPAT 0
+
+#if TENSOR_GGML_COMPAT
+#include "../source/ggml.h"
+#endif
+
+namespace CpuCompute
+{
+	using DirectCompute::TensorShape;
+	using DirectCompute::eDataType;
+
+	__interface iMemoryAllocator
+	{
+		void* allocate( size_t cb, size_t align );
+	};
+	__interface iArenaAllocator : public iMemoryAllocator
+	{
+		void resetArena();
+	};
+
+#if TENSOR_GGML_COMPAT
+	class Tensor;
+	class GgmlTensorView
+	{
+		ggml_tensor tensor;
+	public:
+		GgmlTensorView( const Tensor& t );
+		operator ggml_tensor* ( ) { return &tensor; }
+	};
+#endif
+
+	// A functional equivalent of ggml_tensor structure, designed for use from C++
+	class Tensor : public TensorShape
+	{
+		void* m_data = nullptr;
+
+		eDataType m_type = (eDataType)0xFF;
+
+#if TENSOR_INTERNAL_ALLOC
+		// True when the memory block was allocated internally by this class
+		// In this case, this class does reference counting to support cheap copies.
+		// False when it's owned by someone else, such as iMemoryAllocator object, or a GGML's tensor
+		bool ownsMemory = false;
+		void deallocate();
+#endif
+
+		// Private constructors for fromData() methods
+		Tensor( void* pointer, eDataType type, std::initializer_list<uint32_t> size );
+		Tensor( void* pointer, eDataType type, uint32_t length ) noexcept;
+	public:
+		// Trivial constructors
+		Tensor() = default;
+#if TENSOR_INTERNAL_ALLOC
+		~Tensor()
+		{
+			deallocate();
+		}
+#else
+		~Tensor() = default;
+#endif
+		Tensor( Tensor&& that ) noexcept;
+		void operator=( Tensor&& that ) noexcept;
+		Tensor( const Tensor& that );
+		void operator=( const Tensor& that );
+
+		// Allocate a new tensor
+		HRESULT create( eDataType type, const std::array<uint32_t, 4>& sizeElements, iMemoryAllocator* alloc = nullptr );
+		// Allocate a new tensor
+		HRESULT create( eDataType type, std::initializer_list<uint32_t> sizeElements, iMemoryAllocator* alloc = nullptr );
+		// Attach to pre-existing block of memory, interpreting the data as a dense tensor of the specified type and size
+		HRESULT attach( void* pointer, eDataType type, std::initializer_list<uint32_t> sizeElements );
+		// Attach to pre-existing block of memory, interpret the data as a dense vector of the specified type and length
+		static Tensor fromData( void* pointer, eDataType type, uint32_t length );
+
+		eDataType type() const { return m_type; }
+		void* data() const { return m_data; }
+
+		uint16_t* fp16()
+		{
+			assert( m_type == eDataType::FP16 );
+			assert( nullptr != m_data );
+			return (uint16_t*)m_data;
+		}
+		const uint16_t* fp16() const
+		{
+			assert( m_type == eDataType::FP16 );
+			assert( nullptr != m_data );
+			return (uint16_t*)m_data;
+		}
+		float* fp32()
+		{
+			assert( m_type == eDataType::FP32 );
+			assert( nullptr != m_data );
+			return (float*)m_data;
+		}
+		const float* fp32() const
+		{
+			assert( m_type == eDataType::FP32 );
+			assert( nullptr != m_data );
+			return (float*)m_data;
+		}
+
+		Tensor reshape3d( uint32_t ne0, uint32_t ne1, uint32_t ne2 ) const;
+
+		void setType( eDataType dt )
+		{
+			m_type = dt;
+		}
+		void setDataPointer( void* pv )
+		{
+			m_data = pv;
+		}
+
+#if TENSOR_GGML_COMPAT
+		// Compatibility with GGML's tensors, for testing and lulz
+		Tensor( const ggml_tensor* ggml );
+		ggml_tensor ggml() const;
+
+		operator GgmlTensorView() const
+		{
+			return GgmlTensorView( *this );
+		}
+#endif
+	};
+
+	// A pair of tensors containing weights and biases; apparently, both tensors are of the same shape
+	struct TensorPair
+	{
+		Tensor w, b;
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/TensorCpu.cpp b/Whisper/CPU/TensorCpu.cpp
new file mode 100644
index 0000000..dc34464
--- /dev/null
+++ b/Whisper/CPU/TensorCpu.cpp
@@ -0,0 +1,401 @@
+#include <stdafx.h>
+#include <atomic>
+#include "Tensor.h"
+using namespace CpuCompute;
+
+#if TENSOR_INTERNAL_ALLOC
+namespace
+{
+	// This structure is immediately before the payload of every tensor which has an internally-allocated memory buffer
+	class alignas( 32 ) sTensorMemoryHeader
+	{
+		std::atomic_ptrdiff_t refCounter;
+	public:
+		// Reset the counter to the specified value
+		void reset( ptrdiff_t rc )
+		{
+			refCounter = rc;
+		}
+		// Increment the ref.counter
+		void increment()
+		{
+			refCounter++;
+		}
+		// Decrement the ref.counter, and return true if it reached zero as the result
+		bool decrement()
+		{
+			ptrdiff_t val = --refCounter;
+			assert( val >= 0 );
+			return 0 == val;
+		}
+	};
+
+	inline sTensorMemoryHeader* getMemBlockHeader( void* pv )
+	{
+		assert( nullptr != pv );
+		uint8_t* pb = (uint8_t*)pv;
+		static_assert( sizeof( sTensorMemoryHeader ) == 32 );
+		return (sTensorMemoryHeader*)( pb - sizeof( sTensorMemoryHeader ) );
+	}
+
+	inline void releaseBlock( sTensorMemoryHeader* pointer )
+	{
+		assert( nullptr != pointer );
+		_aligned_free( pointer );
+	}
+
+	inline void* allocateBlock( size_t cb, ptrdiff_t initialRefCounter = 1 )
+	{
+		cb += sizeof( sTensorMemoryHeader );
+		void* pv = _aligned_malloc( cb, 32 );
+		if( nullptr == pv )
+			return nullptr;
+
+		sTensorMemoryHeader* header = (sTensorMemoryHeader*)pv;
+		header->reset( initialRefCounter );
+		return ( (uint8_t*)pv ) + sizeof( sTensorMemoryHeader );
+	}
+}
+
+void Tensor::deallocate()
+{
+	if( ownsMemory && nullptr != m_data )
+	{
+		sTensorMemoryHeader* const header = getMemBlockHeader( m_data );
+		if( header->decrement() )
+		{
+			// This tensor is the last one which had a reference to that block of memory
+			// Release the memory back to the heap
+			releaseBlock( header );
+		}
+	}
+	ownsMemory = false;
+
+	TensorShape::setZero();
+	m_data = nullptr;
+	m_type = (eDataType)0xFF;
+}
+#endif
+
+Tensor::Tensor( const Tensor& that )
+{
+	store( ne, that.sizeVec() );
+	store( nb, that.stridesVec() );
+	m_data = that.m_data;
+	m_type = that.m_type;
+#if TENSOR_INTERNAL_ALLOC
+	if( that.ownsMemory && nullptr != m_data )
+	{
+		getMemBlockHeader( m_data )->increment();
+		ownsMemory = true;
+	}
+	else
+		ownsMemory = false;
+#endif
+}
+
+Tensor::Tensor( Tensor&& that ) noexcept
+{
+	store( ne, that.sizeVec() );
+	store( nb, that.stridesVec() );
+	m_data = that.m_data;
+	m_type = that.m_type;
+#if TENSOR_INTERNAL_ALLOC
+	ownsMemory = that.ownsMemory;
+	that.ownsMemory = false;
+#endif
+	that.m_data = nullptr;
+}
+
+void Tensor::operator=( const Tensor& that )
+{
+	assert( this != &that );
+#if TENSOR_INTERNAL_ALLOC
+	deallocate();
+#endif
+
+	store( ne, that.sizeVec() );
+	store( nb, that.stridesVec() );
+	m_data = that.m_data;
+	m_type = that.m_type;
+#if TENSOR_INTERNAL_ALLOC
+	if( that.ownsMemory && nullptr != m_data )
+	{
+		getMemBlockHeader( m_data )->increment();
+		ownsMemory = true;
+	}
+	else
+		ownsMemory = false;
+#endif
+}
+
+void Tensor::operator=( Tensor&& that ) noexcept
+{
+	assert( this != &that );
+#if TENSOR_INTERNAL_ALLOC
+	deallocate();
+#endif
+	store( ne, that.sizeVec() );
+	store( nb, that.stridesVec() );
+	m_data = that.m_data;
+	m_type = that.m_type;
+	that.m_data = nullptr;
+#if TENSOR_INTERNAL_ALLOC
+	ownsMemory = that.ownsMemory;
+	that.ownsMemory = false;
+#endif
+}
+
+HRESULT Tensor::create( eDataType type, const std::array<uint32_t, 4>& sizeElements, iMemoryAllocator* alloc )
+{
+	const size_t len = (size_t)sizeElements[ 0 ] * sizeElements[ 1 ] * sizeElements[ 2 ] * sizeElements[ 3 ];
+	const size_t cbElement = DirectCompute::elementSize( type );
+	const size_t cb = len * cbElement;
+
+#if TENSOR_INTERNAL_ALLOC
+	deallocate();
+#endif
+
+	store( ne, load( sizeElements ) );
+	TensorShape::setDenseStrides();
+	this->m_type = type;
+
+	if( nullptr != alloc )
+	{
+#if TENSOR_INTERNAL_ALLOC
+		ownsMemory = false;
+#endif
+		m_data = alloc->allocate( cb, 32 );
+		if( nullptr == m_data )
+			return E_OUTOFMEMORY;
+		return S_OK;
+	}
+	else
+	{
+#if TENSOR_INTERNAL_ALLOC
+		m_data = allocateBlock( cb, 1 );
+		if( nullptr == m_data )
+			return E_OUTOFMEMORY;
+		ownsMemory = true;
+		return S_OK;
+#else
+		return E_POINTER;
+#endif
+	}
+}
+
+namespace
+{
+	static HRESULT arrayFromList( std::array<uint32_t, 4>& arr, std::initializer_list<uint32_t> list )
+	{
+		const size_t dims = list.size();
+		if( dims == 0 || dims > 4 )
+			return E_INVALIDARG;
+
+		for( size_t i = 0; i < dims; i++ )
+		{
+			uint32_t u = list.begin()[ i ];
+			if( u == 0 )
+				return E_INVALIDARG;
+			arr[ i ] = u;
+		}
+
+		for( size_t i = dims; i < 4; i++ )
+			arr[ i ] = 1;
+
+		return S_OK;
+	}
+}
+
+HRESULT Tensor::create( eDataType type, std::initializer_list<uint32_t> sizeElements, iMemoryAllocator* alloc )
+{
+	std::array<uint32_t, 4> arr;
+	CHECK( arrayFromList( arr, sizeElements ) );
+
+	return create( type, arr, alloc );
+}
+
+Tensor::Tensor( void* pointer, eDataType type, std::initializer_list<uint32_t> size )
+{
+	if( nullptr == pointer )
+		throw E_POINTER;
+	check( arrayFromList( ne, size ) );
+	TensorShape::setDenseStrides();
+	m_data = pointer;
+	m_type = type;
+#if TENSOR_INTERNAL_ALLOC
+	ownsMemory = false;
+#endif
+}
+
+Tensor::Tensor( void* pointer, eDataType type, uint32_t length ) noexcept
+{
+	// size = [ length, 1, 1, 1 ]
+	const __m128i one = _mm_set1_epi32( 1 );
+	__m128i v = _mm_insert_epi32( one, (int)length, 0 );
+	store( ne, v );
+	// stride = [ 1, length, length, length ]
+	v = _mm_shuffle_epi32( v, _MM_SHUFFLE( 0, 0, 0, 1 ) );
+	store( nb, v );
+
+	m_data = pointer;
+	m_type = type;
+#if TENSOR_INTERNAL_ALLOC
+	ownsMemory = false;
+#endif
+}
+
+Tensor Tensor::fromData( void* pointer, eDataType type, uint32_t length )
+{
+	HRESULT hr = E_UNEXPECTED;
+	if( nullptr != pointer )
+	{
+		if( 0 != length )
+			return Tensor{ pointer, type, length };
+		else
+			hr = E_INVALIDARG;
+	}
+	else
+		hr = E_POINTER;
+	throw hr;
+}
+
+HRESULT Tensor::attach( void* pointer, eDataType type, std::initializer_list<uint32_t> sizeElements )
+{
+	if( nullptr == pointer )
+		return E_POINTER;
+
+	std::array<uint32_t, 4> arr;
+	CHECK( arrayFromList( arr, sizeElements ) );
+
+#if TENSOR_INTERNAL_ALLOC
+	deallocate();
+#endif
+	store( ne, load( arr ) );
+	TensorShape::setDenseStrides();
+
+	m_data = pointer;
+	this->m_type = type;
+#if TENSOR_INTERNAL_ALLOC
+	ownsMemory = false;
+#endif
+	return S_OK;
+}
+
+Tensor Tensor::reshape3d( uint32_t ne0, uint32_t ne1, uint32_t ne2 ) const
+{
+	if( !isContinuous() )
+		throw E_NOTIMPL;
+	if( countElements() != ne0 * ne1 * ne2 )
+		throw E_INVALIDARG;
+
+	Tensor res = *this;
+	res.ne = { ne0, ne1, ne2, 1 };
+	res.setDenseStrides();
+	return res;
+}
+
+#if TENSOR_GGML_COMPAT
+static const __m128i s_maskAlignment16 = _mm_set1_epi64x( 1 );
+static const __m128i s_maskAlignment32 = _mm_set1_epi64x( 3 );
+
+bool isAlignedProperly( __m128i r0, __m128i r1, __m128i mask )
+{
+	__m128i test = _mm_or_si128( r0, r1 );
+	return (bool)_mm_testz_si128( test, mask );
+}
+
+Tensor::Tensor( const ggml_tensor* ggml )
+{
+	store( ne, load16( ggml->ne ) );
+
+	__m128i r0 = load16( (const int*)&ggml->nb[ 0 ] );
+	__m128i r1 = load16( (const int*)&ggml->nb[ 2 ] );
+	// Divide from bytes into elements by right-shifting the 64-bit integers in these vectors
+	switch( ggml->type )
+	{
+	case GGML_TYPE_F16:
+		assert( isAlignedProperly( r0, r1, s_maskAlignment16 ) );
+		r0 = _mm_srli_epi64( r0, 1 );
+		r1 = _mm_srli_epi64( r1, 1 );
+		m_type = eDataType::FP16;
+		break;
+
+	case GGML_TYPE_F32:
+		assert( isAlignedProperly( r0, r1, s_maskAlignment32 ) );
+		r0 = _mm_srli_epi64( r0, 2 );
+		r1 = _mm_srli_epi64( r1, 2 );
+		m_type = eDataType::FP32;
+		break;
+
+	case GGML_TYPE_I32:
+		assert( isAlignedProperly( r0, r1, s_maskAlignment32 ) );
+		r0 = _mm_srli_epi64( r0, 2 );
+		r1 = _mm_srli_epi64( r1, 2 );
+		m_type = eDataType::U32;
+		break;
+
+	default:
+		throw E_INVALIDARG;
+	}
+	// downcast uint64_t into uint32_t in a single vector
+	r0 = _mm_shuffle_epi32( r0, _MM_SHUFFLE( 3, 3, 2, 0 ) );
+	r1 = _mm_shuffle_epi32( r1, _MM_SHUFFLE( 2, 0, 3, 3 ) );
+	store( nb, _mm_blend_epi16( r0, r1, 0b11110000 ) );
+
+	m_data = ggml->data;
+}
+
+ggml_tensor Tensor::ggml() const
+{
+	ggml_tensor res;
+	memset( &res, 0, sizeof( ggml_tensor ) );
+
+	const __m128i size = sizeVec();
+	store16( res.ne, size );
+
+	const __m128i one = _mm_set1_epi32( 1 );
+	const uint32_t maskOnes = (uint32_t)_mm_movemask_ps( _mm_castsi128_ps( _mm_cmpeq_epi32( size, one ) ) );
+	const uint32_t maskNotOnes = maskOnes ^ 0b1111;
+	unsigned long idx;
+	if( _BitScanReverse( &idx, maskNotOnes ) )
+		res.n_dims = (int)idx + 1;
+	else
+		res.n_dims = 0;
+
+	const __m128i strides = stridesVec();
+	// Upcast strides from u32 to u64
+	const __m128i zero = _mm_setzero_si128();
+	__m128i r0 = _mm_unpacklo_epi32( strides, zero );
+	__m128i r1 = _mm_unpackhi_epi32( strides, zero );
+	// Scale from elements into bytes with left shift vector instructions
+	switch( m_type )
+	{
+	case eDataType::FP16:
+		r0 = _mm_slli_epi64( r0, 1 );
+		r1 = _mm_slli_epi64( r1, 1 );
+		res.type = GGML_TYPE_F16;
+		break;
+	case eDataType::FP32:
+		r0 = _mm_slli_epi64( r0, 2 );
+		r1 = _mm_slli_epi64( r1, 2 );
+		res.type = GGML_TYPE_F32;
+		break;
+	case eDataType::U32:
+		r0 = _mm_slli_epi64( r0, 2 );
+		r1 = _mm_slli_epi64( r1, 2 );
+		res.type = GGML_TYPE_I32;
+		break;
+	default:
+		throw OLE_E_BLANK;
+	}
+
+	store16( &res.nb[ 0 ], r0 );
+	store16( &res.nb[ 2 ], r1 );
+
+	res.data = m_data;
+	return res;
+}
+
+GgmlTensorView::GgmlTensorView( const Tensor& t ) : tensor( t.ggml() ) {}
+#endif
\ No newline at end of file
diff --git a/Whisper/CPU/mulMat.cpp b/Whisper/CPU/mulMat.cpp
new file mode 100644
index 0000000..1b6ed25
--- /dev/null
+++ b/Whisper/CPU/mulMat.cpp
@@ -0,0 +1,54 @@
+#include "stdafx.h"
+#include "mulMat.h"
+#include "mulMatImpl.h"
+using namespace CpuCompute;
+
+namespace
+{
+	template<uint8_t panelHeightRegs, uint8_t tileWidthFloats>
+	static HRESULT mulMatImpl( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor )
+	{
+		MulMatImpl<panelHeightRegs, tileWidthFloats> impl{ result, a, b, pfor };
+		return impl.run( pfor );
+	}
+}
+
+HRESULT CpuCompute::mulMat( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor )
+{
+	if( a.type() != eDataType::FP16 )
+		return E_NOTIMPL;
+	if( b.type() != eDataType::FP32 )
+		return E_NOTIMPL;
+
+	// return mulMatImpl<1, 1>( result, a, b, pfor );
+
+	if( b.ne[ 1 ] == 1 )
+	{
+		// Multiplying by a single row
+		if( a.ne[ 1 ] >= 32 )
+			return mulMatImpl<4, 1>( result, a, b, pfor );
+		else
+			return mulMatImpl<1, 1>( result, a, b, pfor );
+	}
+	else if( b.ne[ 1 ] == 2 )
+	{
+		if( a.ne[ 1 ] >= 32 )
+			return mulMatImpl<4, 2>( result, a, b, pfor );
+		else
+			return mulMatImpl<1, 2>( result, a, b, pfor );
+	}
+	else if( b.ne[ 1 ] == 3 )
+	{
+		if( a.ne[ 1 ] >= 16 )
+			return mulMatImpl<2, 3>( result, a, b, pfor );
+		else
+			return mulMatImpl<1, 3>( result, a, b, pfor );
+	}
+	else
+	{
+		if( a.ne[ 1 ] >= 16 )
+			return mulMatImpl<2, 4>( result, a, b, pfor );
+		else
+			return mulMatImpl<1, 4>( result, a, b, pfor );
+	}
+}
\ No newline at end of file
diff --git a/Whisper/CPU/mulMat.h b/Whisper/CPU/mulMat.h
new file mode 100644
index 0000000..36e56fe
--- /dev/null
+++ b/Whisper/CPU/mulMat.h
@@ -0,0 +1,17 @@
+#pragma once
+#include "ParallelForRunner.h"
+#include "Tensor.h"
+
+namespace CpuCompute
+{
+	HRESULT mulMat( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor );
+}
+
+#if TENSOR_GGML_COMPAT
+#include "../source/ggml.h"
+inline HRESULT mulMat( ggml_tensor* result, const ggml_tensor* a, const ggml_tensor* b, CpuCompute::ParallelForRunner& pfor )
+{
+	CpuCompute::Tensor r{ result }, lhs{ a }, rhs{ b };
+	return CpuCompute::mulMat( r, lhs, rhs, pfor );
+}
+#endif
\ No newline at end of file
diff --git a/Whisper/CPU/mulMat.kernel.hpp b/Whisper/CPU/mulMat.kernel.hpp
new file mode 100644
index 0000000..80b51a0
--- /dev/null
+++ b/Whisper/CPU/mulMat.kernel.hpp
@@ -0,0 +1,742 @@
+#pragma once
+#include <stdint.h>
+#include <immintrin.h>
+#include "simdUtils.h"
+
+template<uint8_t panelHeightRegs, uint8_t tileWidthFloats>
+struct ResultTile
+{
+	static constexpr size_t totalRegs = (size_t)(tileWidthFloats)*panelHeightRegs;
+	std::array<__m256, totalRegs> arr;
+
+	template<size_t idx>
+	__forceinline void fmadd( __m256 a, __m256 b )
+	{
+		arr[ idx ] = _mm256_fmadd_ps( a, b, arr[ idx ] );
+	}
+	__forceinline void kernel( const std::array<__m256, panelHeightRegs>& panel, const float* rsi, size_t stride );
+	__forceinline void kernelPartial( const std::array<__m256, panelHeightRegs>& panel, const float* rsi, size_t stride, size_t rem )
+	{
+		throw E_UNEXPECTED;
+	}
+	__forceinline void store( float* rdi, size_t w, size_t h, size_t stride ) const;
+};
+
+#pragma region setZero functions
+__forceinline void setZero( std::array<__m256, 1>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+}
+__forceinline void setZero( std::array<__m256, 2>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+	dest[ 1 ] = _mm256_setzero_ps();
+}
+__forceinline void setZero( std::array<__m256, 3>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+	dest[ 1 ] = _mm256_setzero_ps();
+	dest[ 2 ] = _mm256_setzero_ps();
+}
+__forceinline void setZero( std::array<__m256, 4>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+	dest[ 1 ] = _mm256_setzero_ps();
+	dest[ 2 ] = _mm256_setzero_ps();
+	dest[ 3 ] = _mm256_setzero_ps();
+}
+__forceinline void setZero( std::array<__m256, 6>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+	dest[ 1 ] = _mm256_setzero_ps();
+	dest[ 2 ] = _mm256_setzero_ps();
+	dest[ 3 ] = _mm256_setzero_ps();
+	dest[ 4 ] = _mm256_setzero_ps();
+	dest[ 5 ] = _mm256_setzero_ps();
+}
+__forceinline void setZero( std::array<__m256, 8>& dest )
+{
+	dest[ 0 ] = _mm256_setzero_ps();
+	dest[ 1 ] = _mm256_setzero_ps();
+	dest[ 2 ] = _mm256_setzero_ps();
+	dest[ 3 ] = _mm256_setzero_ps();
+	dest[ 4 ] = _mm256_setzero_ps();
+	dest[ 5 ] = _mm256_setzero_ps();
+	dest[ 6 ] = _mm256_setzero_ps();
+	dest[ 7 ] = _mm256_setzero_ps();
+}
+#pragma endregion
+
+#pragma region Micro-kernels
+__forceinline void ResultTile<1, 1>::kernel( const std::array<__m256, 1>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+}
+__forceinline void ResultTile<1, 2>::kernel( const std::array<__m256, 1>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<1>( panel[ 0 ], b );
+}
+__forceinline void ResultTile<1, 2>::kernelPartial( const std::array<__m256, 1>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( 1 == rem );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+}
+__forceinline void ResultTile<1, 3>::kernel( const std::array<__m256, 1>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<1>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride * 2 );
+	fmadd<2>( panel[ 0 ], b );
+}
+__forceinline void ResultTile<1, 3>::kernelPartial( const std::array<__m256, 1>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( rem > 0 && rem < 3 );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	if( rem > 1 )
+	{
+		b = _mm256_broadcast_ss( rsi + stride );
+		fmadd<1>( panel[ 0 ], b );
+	}
+}
+
+__forceinline void ResultTile<1, 4>::kernel( const std::array<__m256, 1>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<1>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride * 2 );
+	fmadd<2>( panel[ 0 ], b );
+	b = _mm256_broadcast_ss( rsi + stride * 3 );
+	fmadd<3>( panel[ 0 ], b );
+}
+__forceinline void ResultTile<1, 4>::kernelPartial( const std::array<__m256, 1>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( rem > 0 && rem < 4 );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+
+	switch( rem )
+	{
+	case 3:
+		b = _mm256_broadcast_ss( rsi + stride * 2 );
+		fmadd<2>( panel[ 0 ], b );
+	case 2:
+		b = _mm256_broadcast_ss( rsi + stride );
+		fmadd<1>( panel[ 0 ], b );
+	}
+}
+__forceinline void ResultTile<4, 1>::kernel( const std::array<__m256, 4>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+	fmadd<2>( panel[ 2 ], b );
+	fmadd<3>( panel[ 3 ], b );
+}
+__forceinline void ResultTile<2, 4>::kernel( const std::array<__m256, 2>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<2>( panel[ 0 ], b );
+	fmadd<3>( panel[ 1 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride * 2 );
+	fmadd<4>( panel[ 0 ], b );
+	fmadd<5>( panel[ 1 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride * 3 );
+	fmadd<6>( panel[ 0 ], b );
+	fmadd<7>( panel[ 1 ], b );
+}
+
+__forceinline void ResultTile<2, 4>::kernelPartial( const std::array<__m256, 2>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( rem > 0 && rem < 4 );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+
+	switch( rem )
+	{
+	case 3:
+		b = _mm256_broadcast_ss( rsi + stride * 2 );
+		fmadd<4>( panel[ 0 ], b );
+		fmadd<5>( panel[ 1 ], b );
+	case 2:
+		b = _mm256_broadcast_ss( rsi + stride );
+		fmadd<2>( panel[ 0 ], b );
+		fmadd<3>( panel[ 1 ], b );
+	}
+}
+
+__forceinline void ResultTile<2, 3>::kernel( const std::array<__m256, 2>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<2>( panel[ 0 ], b );
+	fmadd<3>( panel[ 1 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride * 2 );
+	fmadd<4>( panel[ 0 ], b );
+	fmadd<5>( panel[ 1 ], b );
+}
+__forceinline void ResultTile<2, 3>::kernelPartial( const std::array<__m256, 2>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( rem > 0 && rem < 3 );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+	if( rem > 1 )
+	{
+		b = _mm256_broadcast_ss( rsi + stride );
+		fmadd<2>( panel[ 0 ], b );
+		fmadd<3>( panel[ 1 ], b );
+	}
+}
+
+__forceinline void ResultTile<4, 2>::kernel( const std::array<__m256, 4>& panel, const float* rsi, size_t stride )
+{
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+	fmadd<2>( panel[ 2 ], b );
+	fmadd<3>( panel[ 3 ], b );
+
+	b = _mm256_broadcast_ss( rsi + stride );
+	fmadd<4>( panel[ 0 ], b );
+	fmadd<5>( panel[ 1 ], b );
+	fmadd<6>( panel[ 2 ], b );
+	fmadd<7>( panel[ 3 ], b );
+}
+__forceinline void ResultTile<4, 2>::kernelPartial( const std::array<__m256, 4>& panel, const float* rsi, size_t stride, size_t rem )
+{
+	assert( 1 == rem );
+	__m256 b = _mm256_broadcast_ss( rsi );
+	fmadd<0>( panel[ 0 ], b );
+	fmadd<1>( panel[ 1 ], b );
+	fmadd<2>( panel[ 2 ], b );
+	fmadd<3>( panel[ 3 ], b );
+}
+#pragma endregion
+
+#pragma region Loads
+// This function should compile into a single `vcvtph2ps` instruction, with memory operand
+__forceinline __m256 loadUpcasted( const uint16_t* rsi )
+{
+	__m128i i = _mm_load_si128( ( const __m128i* )rsi );
+	return _mm256_cvtph_ps( i );
+}
+
+// We loading the panel from the temporary buffer.
+// For this reason, we don't need to handle remainders, the code which made the buffer wrote zeros into the remainder elements
+// We can even use aligned load instructions.
+__forceinline void loadPanel( const uint16_t* rsi, std::array<__m256, 1>& dest )
+{
+	dest[ 0 ] = loadUpcasted( rsi );
+}
+__forceinline void loadPanel( const uint16_t* rsi, std::array<__m256, 2>& dest )
+{
+	dest[ 0 ] = loadUpcasted( rsi );
+	dest[ 1 ] = loadUpcasted( rsi + 8 );
+}
+__forceinline void loadPanel( const uint16_t* rsi, std::array<__m256, 3>& dest )
+{
+	dest[ 0 ] = loadUpcasted( rsi );
+	dest[ 1 ] = loadUpcasted( rsi + 8 );
+	dest[ 2 ] = loadUpcasted( rsi + 8 * 2 );
+}
+__forceinline void loadPanel( const uint16_t* rsi, std::array<__m256, 4>& dest )
+{
+	dest[ 0 ] = loadUpcasted( rsi );
+	dest[ 1 ] = loadUpcasted( rsi + 8 );
+	dest[ 2 ] = loadUpcasted( rsi + 8 * 2 );
+	dest[ 3 ] = loadUpcasted( rsi + 8 * 3 );
+}
+#pragma endregion
+
+#pragma region Stores
+__forceinline void ResultTile<1, 1>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h == 1 && w > 0 && w <= 8 );
+	if( w == 8 )
+		_mm256_storeu_ps( rdi, arr[ 0 ] );
+	else
+	{
+		const __m256i mask = loadTailMaskInt( w );
+		_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+	}
+}
+
+__forceinline void ResultTile<1, 2>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 2 && w <= 8 );
+	if( w == 8 )
+	{
+		switch( h )
+		{
+		case 2:
+			_mm256_storeu_ps( rdi + stride, arr[ 1 ] );
+		case 1:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+		}
+	}
+	else
+	{
+		const __m256i mask = loadTailMaskInt( w );
+		switch( h )
+		{
+		case 2:
+			_mm256_maskstore_ps( rdi + stride, mask, arr[ 1 ] );
+		case 1:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+		}
+	}
+}
+
+__forceinline void ResultTile<1, 3>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 3 && w <= 8 );
+	if( w == 8 )
+	{
+		switch( h )
+		{
+		case 3:
+			_mm256_storeu_ps( rdi + stride * 2, arr[ 2 ] );
+		case 2:
+			_mm256_storeu_ps( rdi + stride, arr[ 1 ] );
+		case 1:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+		}
+	}
+	else
+	{
+		const __m256i mask = loadTailMaskInt( w );
+		switch( h )
+		{
+		case 3:
+			_mm256_maskstore_ps( rdi + stride * 2, mask, arr[ 2 ] );
+		case 2:
+			_mm256_maskstore_ps( rdi + stride, mask, arr[ 1 ] );
+		case 1:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+		}
+	}
+}
+
+__forceinline void ResultTile<1, 4>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 4 && w <= 8 );
+
+	if( w == 8 )
+	{
+		switch( h )
+		{
+		case 4:
+			_mm256_storeu_ps( rdi + stride * 3, arr[ 3 ] );
+		case 3:
+			_mm256_storeu_ps( rdi + stride * 2, arr[ 2 ] );
+		case 2:
+			_mm256_storeu_ps( rdi + stride, arr[ 1 ] );
+		case 1:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+		}
+	}
+	else
+	{
+		const __m256i mask = loadTailMaskInt( w );
+		switch( h )
+		{
+		case 4:
+			_mm256_maskstore_ps( rdi + stride * 3, mask, arr[ 3 ] );
+		case 3:
+			_mm256_maskstore_ps( rdi + stride * 2, mask, arr[ 2 ] );
+		case 2:
+			_mm256_maskstore_ps( rdi + stride, mask, arr[ 1 ] );
+		case 1:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+		}
+	}
+}
+
+__forceinline void ResultTile<4, 1>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h == 1 && w > 0 && w <= 32 );
+	if( w == 32 )
+	{
+		// 4 complete vectors, this branch is very likely to be taken
+		_mm256_storeu_ps( rdi, arr[ 0 ] );
+		_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+		_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+		_mm256_storeu_ps( rdi + 8 * 3, arr[ 3 ] );
+	}
+	else
+	{
+		const size_t rem = w % 8;
+		const __m256i mask = loadTailMaskInt<false>( rem );
+		const size_t completeVectors = w / 8;
+		const size_t key = ( completeVectors << 1 ) | ( ( 0 == rem ) ? 0 : 1 );
+		switch( key )
+		{
+		case 1:	// 0 complete vectors + remainder
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			break;
+		case 2:	// 1 complete vector
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			break;
+		case 3:	// 1 complete vector + remainder
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			break;
+		case 4:	// 2 complete vectors
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			break;
+		case 5:	// 2 complete vectors + remainder
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_maskstore_ps( rdi + 8 * 2, mask, arr[ 2 ] );
+			break;
+		case 6:	// 3 complete vectors
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+			break;
+		case 7:	// 3 complete vectors + remainder
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi + 8 * 3, mask, arr[ 3 ] );
+			break;
+		default:
+			throw E_UNEXPECTED;
+		}
+	}
+}
+__forceinline void ResultTile<4, 2>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 2 && w <= 32 );
+	const bool twoRows = h == 2;
+	float* const rdi1 = rdi + stride;
+	if( w == 32 )
+	{
+		_mm256_storeu_ps( rdi, arr[ 0 ] );
+		_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+		_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+		_mm256_storeu_ps( rdi + 8 * 3, arr[ 3 ] );
+
+		if( twoRows )
+		{
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 5 ] );
+			_mm256_storeu_ps( rdi1 + 8 * 2, arr[ 6 ] );
+			_mm256_storeu_ps( rdi1 + 8 * 3, arr[ 7 ] );
+		}
+	}
+	else
+	{
+		const size_t rem = w % 8;
+		const __m256i mask = loadTailMaskInt<false>( rem );
+		const size_t completeVectors = w / 8;
+		// Lowest bit: remainder
+		// Next bit: set when storing 2 rows
+		// Next 2 bits: count of complete vectors in X direction, [ 0..3 ]
+		const size_t key = ( completeVectors << 2 ) | ( ( 0 == rem ) ? 0 : 1 ) | ( twoRows ? 2 : 0 );
+		switch( key )
+		{
+		case 1:	// 0 complete vectors + remainder, 1 row
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			break;
+		case 3:	// 0 complete vectors + remainder, 2 rows
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 4 ] );
+			break;
+		case 4:	// 1 complete vector, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			break;
+		case 5:	// 1 complete vector + remainder, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			break;
+		case 6:	// 1 complete vector, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			break;
+		case 7:	// 1 complete vector + remainder, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 5 ] );
+			break;
+		case 8:	// 2 complete vectors, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			break;
+		case 9:	// 2 complete vectors + remainder, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_maskstore_ps( rdi + 8 * 2, mask, arr[ 2 ] );
+			break;
+		case 10:	// 2 complete vectors, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 5 ] );
+			break;
+		case 11:	// 2 complete vectors + remainder, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_maskstore_ps( rdi + 8 * 2, mask, arr[ 2 ] );
+
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 5 ] );
+			_mm256_maskstore_ps( rdi1 + 8 * 2, mask, arr[ 6 ] );
+			break;
+		case 12:	// 3 complete vectors, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+			break;
+		case 13:	// 3 complete vectors + remainder, 1 row
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi + 8 * 3, mask, arr[ 3 ] );
+			break;
+		case 14:	// 3 complete vectors, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 5 ] );
+			_mm256_storeu_ps( rdi1 + 8 * 2, arr[ 6 ] );
+			break;
+		case 15:	// 3 complete vectors + remainder, 2 rows
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+			_mm256_storeu_ps( rdi + 8 * 2, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi + 8 * 3, mask, arr[ 3 ] );
+
+			_mm256_storeu_ps( rdi1, arr[ 4 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 5 ] );
+			_mm256_storeu_ps( rdi1 + 8 * 2, arr[ 6 ] );
+			_mm256_maskstore_ps( rdi1 + 8 * 3, mask, arr[ 7 ] );
+			break;
+		default:
+			throw E_UNEXPECTED;
+		}
+	}
+}
+
+__forceinline void ResultTile<2, 4>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 4 && w <= 16 );
+	h--;
+	float* const rdi1 = rdi + stride;
+	float* const rdi2 = rdi + stride * 2;
+	float* const rdi3 = rdi + stride * 3;
+
+	if( w == 16 )
+	{
+		switch( h )
+		{
+		case 3:
+			_mm256_storeu_ps( rdi3, arr[ 6 ] );
+			_mm256_storeu_ps( rdi3 + 8, arr[ 7 ] );
+		case 2:
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_storeu_ps( rdi2 + 8, arr[ 5 ] );
+		case 1:
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 3 ] );
+		case 0:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+		}
+	}
+	else
+	{
+		const size_t rem = w % 8;
+		const __m256i mask = loadTailMaskInt<false>( rem );
+		// 0 for partial first vector, 1 for exactly 1 complete vector, 2 for 1 complete vector with remainder
+		const size_t partialCase = ( w < 8 ) ? 0 : ( ( w == 8 ) ? 1 : 2 );
+		// Merge into a single integer for the switch statement
+		const size_t key = partialCase + h * 3;
+
+		switch( key )
+		{
+			// h = 1
+		case 0:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			break;
+		case 1:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			break;
+		case 2:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			break;
+			// h = 2
+		case 3:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 2 ] );
+			break;
+		case 4:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			break;
+		case 5:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 3 ] );
+			break;
+			// h = 3
+		case 6:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi2, mask, arr[ 4 ] );
+			break;
+		case 7:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			break;
+		case 8:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 3 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_maskstore_ps( rdi2 + 8, mask, arr[ 5 ] );
+			break;
+			// h = 4
+		case 9:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi2, mask, arr[ 4 ] );
+			_mm256_maskstore_ps( rdi3, mask, arr[ 6 ] );
+			break;
+		case 10:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_storeu_ps( rdi3, arr[ 6 ] );
+			break;
+		case 11:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 3 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_maskstore_ps( rdi2 + 8, mask, arr[ 5 ] );
+			_mm256_storeu_ps( rdi3, arr[ 6 ] );
+			_mm256_maskstore_ps( rdi3 + 8, mask, arr[ 7 ] );
+			break;
+		default:
+			throw E_UNEXPECTED;
+		}
+	}
+}
+
+__forceinline void ResultTile<2, 3>::store( float* rdi, size_t w, size_t h, size_t stride ) const
+{
+	assert( h > 0 && w > 0 && h <= 3 && w <= 16 );
+	float* const rdi1 = rdi + stride;
+	float* const rdi2 = rdi + stride * 2;
+	h--;
+
+	if( w == 16 )
+	{
+		switch( h )
+		{
+		case 2:
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_storeu_ps( rdi2 + 8, arr[ 5 ] );
+		case 1:
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_storeu_ps( rdi1 + 8, arr[ 3 ] );
+		case 0:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi + 8, arr[ 1 ] );
+		}
+	}
+	else
+	{
+		const size_t rem = w % 8;
+		const __m256i mask = loadTailMaskInt<false>( rem );
+		// 0 for partial first vector, 1 for exactly 1 complete vector, 2 for 1 complete vector with remainder
+		const size_t partialCase = ( w < 8 ) ? 0 : ( ( w == 8 ) ? 1 : 2 );
+		// Merge into a single integer for the switch statement
+		const size_t key = partialCase + h * 3;
+
+		switch( key )
+		{
+			// h = 1
+		case 0:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			break;
+		case 1:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			break;
+		case 2:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			break;
+			// h = 2
+		case 3:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 2 ] );
+			break;
+		case 4:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			break;
+		case 5:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 3 ] );
+			break;
+			// h = 3
+		case 6:
+			_mm256_maskstore_ps( rdi, mask, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi1, mask, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi2, mask, arr[ 4 ] );
+			break;
+		case 7:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			break;
+		case 8:
+			_mm256_storeu_ps( rdi, arr[ 0 ] );
+			_mm256_maskstore_ps( rdi + 8, mask, arr[ 1 ] );
+			_mm256_storeu_ps( rdi1, arr[ 2 ] );
+			_mm256_maskstore_ps( rdi1 + 8, mask, arr[ 3 ] );
+			_mm256_storeu_ps( rdi2, arr[ 4 ] );
+			_mm256_maskstore_ps( rdi2 + 8, mask, arr[ 5 ] );
+			break;
+		default:
+			throw E_UNEXPECTED;
+		}
+	}
+}
+#pragma endregion
\ No newline at end of file
diff --git a/Whisper/CPU/mulMatImpl.avx2.cpp b/Whisper/CPU/mulMatImpl.avx2.cpp
new file mode 100644
index 0000000..b15ae63
--- /dev/null
+++ b/Whisper/CPU/mulMatImpl.avx2.cpp
@@ -0,0 +1,362 @@
+#include "stdafx.h"
+#include "mulMatImpl.h"
+#include <immintrin.h>
+#include "mulMatUtils.hpp"
+using namespace CpuCompute;
+
+namespace
+{
+	constexpr size_t prefetchBytes = 96;
+	constexpr int prefetchHint = _MM_HINT_T0;
+
+	constexpr size_t maskAlign16 = ~(size_t)15;
+
+	__forceinline __m256i load( const void* rsi )
+	{
+		return _mm256_loadu_si256( ( const __m256i* )rsi );
+	}
+
+#define TRANSPOSE_8X16()                           \
+                                                   \
+	__m256i t0 = _mm256_unpacklo_epi16( r0, r1 );  \
+	__m256i t1 = _mm256_unpackhi_epi16( r0, r1 );  \
+	__m256i t2 = _mm256_unpacklo_epi16( r2, r3 );  \
+	__m256i t3 = _mm256_unpackhi_epi16( r2, r3 );  \
+	__m256i t4 = _mm256_unpacklo_epi16( r4, r5 );  \
+	__m256i t5 = _mm256_unpackhi_epi16( r4, r5 );  \
+	__m256i t6 = _mm256_unpacklo_epi16( r6, r7 );  \
+	__m256i t7 = _mm256_unpackhi_epi16( r6, r7 );  \
+                                                   \
+	r0 = _mm256_unpacklo_epi32( t0, t2 );          \
+	r1 = _mm256_unpackhi_epi32( t0, t2 );          \
+	r2 = _mm256_unpacklo_epi32( t1, t3 );          \
+	r3 = _mm256_unpackhi_epi32( t1, t3 );          \
+	r4 = _mm256_unpacklo_epi32( t4, t6 );          \
+	r5 = _mm256_unpackhi_epi32( t4, t6 );          \
+	r6 = _mm256_unpacklo_epi32( t5, t7 );          \
+	r7 = _mm256_unpackhi_epi32( t5, t7 );          \
+                                                   \
+	t0 = _mm256_unpacklo_epi64( r0, r4 );          \
+	t1 = _mm256_unpackhi_epi64( r0, r4 );          \
+	t2 = _mm256_unpacklo_epi64( r1, r5 );          \
+	t3 = _mm256_unpackhi_epi64( r1, r5 );          \
+	t4 = _mm256_unpacklo_epi64( r2, r6 );          \
+	t5 = _mm256_unpackhi_epi64( r2, r6 );          \
+	t6 = _mm256_unpacklo_epi64( r3, r7 );          \
+	t7 = _mm256_unpackhi_epi64( r3, r7 )
+
+	__forceinline void storeLow( void* rdi, __m256i v )
+	{
+		__m128i i = _mm256_castsi256_si128( v );
+		_mm_store_si128( ( __m128i* )rdi, i );
+	}
+
+#define STORE_8X16_LOW()                     \
+	storeLow( rdi, t0 );                     \
+	storeLow( rdi + destStride, t1 );        \
+	storeLow( rdi + destStride * 2, t2 );    \
+	rdi += destStride * 8;                   \
+	storeLow( rdiMid, t3 );                  \
+	storeLow( rdiMid + destStride, t4 );     \
+	storeLow( rdiMid + destStride * 2, t5 ); \
+	rdiMid += destStride * 8;                \
+	storeLow( rdiLast, t6 );                 \
+	storeLow( rdiLast + destStride, t7 );    \
+	rdiLast += destStride * 8
+
+	__forceinline void storeHigh( void* rdi, __m256i v )
+	{
+		__m128i i = _mm256_extracti128_si256( v, 1 );
+		_mm_store_si128( ( __m128i* )rdi, i );
+	}
+
+#define STORE_8X16_HIGH()                     \
+	storeHigh( rdi, t0 );                     \
+	storeHigh( rdi + destStride, t1 );        \
+	storeHigh( rdi + destStride * 2, t2 );    \
+	rdi += destStride * 8;                    \
+	storeHigh( rdiMid, t3 );                  \
+	storeHigh( rdiMid + destStride, t4 );     \
+	storeHigh( rdiMid + destStride * 2, t5 ); \
+	rdiMid += destStride * 8;                 \
+	storeHigh( rdiLast, t6 );                 \
+	storeHigh( rdiLast + destStride, t7 );    \
+	rdiLast += destStride * 8
+
+	__forceinline void prefetch( const uint8_t* p )
+	{
+		_mm_prefetch( (const char*)p, prefetchHint );
+	}
+
+	__forceinline void transpose8Avx2( uint16_t* rdiWords, size_t w, const uint16_t* rsiWords, size_t sourceStride, size_t destStride )
+	{
+		assert( 0 == ( (size_t)rdiWords ) % 16 );
+		assert( 0 == destStride % 8 );
+		assert( w <= sourceStride );
+
+		// Scale strides to bytes, and cast the pointers
+		sourceStride *= 2;
+		destStride *= 2;
+		uint8_t* rdi = (uint8_t*)rdiWords;
+		const uint8_t* rsi = (const uint8_t*)rsiWords;
+
+		const uint8_t* const rsiEndAligned = rsi + ( w & maskAlign16 ) * 2;
+		const uint8_t* const rsiEnd = rsi + w * 2;
+		const uint8_t* rsiMid = rsi + sourceStride * 3;
+		const uint8_t* rsiLast = rsi + sourceStride * 6;
+		uint8_t* rdiMid = rdi + destStride * 3;
+		uint8_t* rdiLast = rdi + destStride * 6;
+
+		while( rsi < rsiEndAligned )
+		{
+			// Load 16x8 block into 8 registers
+			__m256i r0 = load( rsi );
+			__m256i r1 = load( rsi + sourceStride );
+			__m256i r2 = load( rsi + sourceStride * 2 );
+			rsi += 32;
+			__m256i r3 = load( rsiMid );
+			__m256i r4 = load( rsiMid + sourceStride );
+			__m256i r5 = load( rsiMid + sourceStride * 2 );
+			rsiMid += 32;
+			__m256i r6 = load( rsiLast );
+			__m256i r7 = load( rsiLast + sourceStride );
+			rsiLast += 32;
+
+			// Transpose FP16 values in registers
+			TRANSPOSE_8X16();
+
+			// Store
+			STORE_8X16_LOW();
+			STORE_8X16_HIGH();
+
+			if constexpr( prefetchBytes > 0 )
+			{
+				if( rsi + prefetchBytes < rsiEnd )
+				{
+					prefetch( rsi + prefetchBytes );
+					prefetch( rsi + sourceStride + prefetchBytes );
+					prefetch( rsi + sourceStride * 2 + prefetchBytes );
+					prefetch( rsiMid + prefetchBytes );
+					prefetch( rsiMid + sourceStride + prefetchBytes );
+					prefetch( rsiMid + sourceStride * 2 + prefetchBytes );
+					prefetch( rsiLast + prefetchBytes );
+					prefetch( rsiLast + sourceStride + prefetchBytes );
+				}
+			}
+		}
+
+		if( rsi < rsiEnd )
+		{
+			// Loading 8 elements into corresponding lanes of 8 vectors
+			// This way there's no data dependencies between these load instructions
+			// Out of order execution should hopefully do it's magic in the CPU, running all these loads in parallel.
+			__m128i r0;
+			__m128i r1 = _mm_setzero_si128();
+			__m128i r2 = _mm_setzero_si128();
+			__m128i r3 = _mm_setzero_si128();
+			__m128i r4 = _mm_setzero_si128();
+			__m128i r5 = _mm_setzero_si128();
+			__m128i r6 = _mm_setzero_si128();
+			__m128i r7 = _mm_setzero_si128();
+
+			__m128i t0, t1, t2, t3, t4, t5, t6;
+
+#pragma loop( no_vector )
+			while( rsi < rsiEnd )
+			{
+				r0 = _mm_cvtsi32_si128( *(const uint16_t*)rsi );
+				r1 = _mm_insert_epi16( r1, *(const int16_t*)( rsi + sourceStride ), 1 );
+				r2 = _mm_insert_epi16( r2, *(const int16_t*)( rsi + sourceStride * 2 ), 2 );
+				rsi += 2;
+				r3 = _mm_insert_epi16( r3, *(const int16_t*)( rsiMid ), 3 );
+				r4 = _mm_insert_epi16( r4, *(const int16_t*)( rsiMid + sourceStride ), 4 );
+				r5 = _mm_insert_epi16( r5, *(const int16_t*)( rsiMid + sourceStride * 2 ), 5 );
+				rsiMid += 2;
+				r6 = _mm_insert_epi16( r6, *(const int16_t*)( rsiLast ), 6 );
+				r7 = _mm_insert_epi16( r7, *(const int16_t*)( rsiLast + sourceStride ), 7 );
+				rsiLast += 2;
+
+				// Bitwise operations are pretty fast, AMD Zen3 CPU can run 4 of them every clock cycle
+				// Combine 8 vectors into one
+				t0 = _mm_or_si128( r0, r1 );
+				t1 = _mm_or_si128( r2, r3 );
+				t2 = _mm_or_si128( r4, r5 );
+				t3 = _mm_or_si128( r6, r7 );
+
+				t4 = _mm_or_si128( t0, t1 );
+				t5 = _mm_or_si128( t2, t3 );
+
+				t6 = _mm_or_si128( t4, t5 );
+				// Store 8 FP16 values, the destination is aligned
+				_mm_store_si128( ( __m128i* )rdi, t6 );
+				rdi += destStride;
+			}
+		}
+	}
+
+	__forceinline void transpose8PartialAvx2( uint16_t* rdiWords, size_t w, size_t h, const uint16_t* rsiWords, size_t sourceStride, size_t destStride )
+	{
+		assert( 0 == ( (size_t)rdiWords ) % 16 );
+		assert( 0 == destStride % 8 );
+		assert( w <= sourceStride );
+		assert( h > 0 && h < 8 );
+
+		// Scale strides to bytes, and cast the pointers
+		sourceStride *= 2;
+		destStride *= 2;
+		uint8_t* rdi = (uint8_t*)rdiWords;
+		const uint8_t* rsi = (const uint8_t*)rsiWords;
+
+		const uint8_t* const rsiEndAligned = rsi + ( w & maskAlign16 ) * 2;
+		const uint8_t* const rsiEnd = rsi + w * 2;
+		const uint8_t* rsiMid = rsi + sourceStride * 3;
+		const uint8_t* rsiLast = rsi + sourceStride * 6;
+		uint8_t* rdiMid = rdi + destStride * 3;
+		uint8_t* rdiLast = rdi + destStride * 6;
+
+		while( rsi < rsiEndAligned )
+		{
+			// Load the block into 8 registers, set unused rows to zero
+			__m256i r0 = load( rsi );
+			__m256i r1 = _mm256_setzero_si256();
+			__m256i r2 = _mm256_setzero_si256();
+			__m256i r3 = _mm256_setzero_si256();
+			__m256i r4 = _mm256_setzero_si256();
+			__m256i r5 = _mm256_setzero_si256();
+			__m256i r6 = _mm256_setzero_si256();
+			// These branches, whether direct or indirect, are very predictable: same outcome for all iterations of the outer loop
+			switch( h )
+			{
+			case 7:
+				r6 = load( rsiLast );
+			case 6:
+				r5 = load( rsiMid + sourceStride * 2 );
+			case 5:
+				r4 = load( rsiMid + sourceStride );
+			case 4:
+				r3 = load( rsiMid );
+			case 3:
+				r2 = load( rsi + sourceStride * 2 );
+			case 2:
+				r1 = load( rsi + sourceStride );
+			}
+			rsi += 32;
+			rsiMid += 32;
+			rsiLast += 32;
+
+			__m256i r7 = _mm256_setzero_si256();
+
+			// Transpose FP16 values in registers
+			TRANSPOSE_8X16();
+
+			// Store
+			STORE_8X16_LOW();
+
+			STORE_8X16_HIGH();
+		}
+
+		if( rsi < rsiEnd )
+		{
+			// Loading 8 elements into corresponding lanes of 8 vectors
+			// This way there's no data dependencies between these load instructions
+			// Out of order execution should hopefully do it's magic in the CPU, running all these loads in parallel.
+			__m128i r0;
+			__m128i r1 = _mm_setzero_si128();
+			__m128i r2 = _mm_setzero_si128();
+			__m128i r3 = _mm_setzero_si128();
+			__m128i r4 = _mm_setzero_si128();
+			__m128i r5 = _mm_setzero_si128();
+			__m128i r6 = _mm_setzero_si128();
+
+			__m128i t0, t1, t2, t3, t4, t5;
+
+#pragma loop( no_vector )
+			while( rsi < rsiEnd )
+			{
+				r0 = _mm_cvtsi32_si128( *(const uint16_t*)rsi );
+
+				switch( h )
+				{
+				case 7:
+					r6 = _mm_insert_epi16( r6, *(const int16_t*)( rsiLast ), 6 );
+				case 6:
+					r5 = _mm_insert_epi16( r5, *(const int16_t*)( rsiMid + sourceStride * 2 ), 5 );
+				case 5:
+					r4 = _mm_insert_epi16( r4, *(const int16_t*)( rsiMid + sourceStride ), 4 );
+				case 4:
+					r3 = _mm_insert_epi16( r3, *(const int16_t*)( rsiMid ), 3 );
+				case 3:
+					r2 = _mm_insert_epi16( r2, *(const int16_t*)( rsi + sourceStride * 2 ), 2 );
+				case 2:
+					r1 = _mm_insert_epi16( r1, *(const int16_t*)( rsi + sourceStride ), 1 );
+				}
+				rsi += 2;
+				rsiMid += 2;
+				rsiLast += 2;
+
+				// Bitwise operations are pretty fast, AMD Zen3 CPU can run 4 of them every clock cycle
+				// Combine 7 vectors into one
+				t0 = _mm_or_si128( r0, r1 );
+				t1 = _mm_or_si128( r2, r3 );
+				t2 = _mm_or_si128( r4, r5 );
+
+				t3 = _mm_or_si128( t0, t1 );
+				t4 = _mm_or_si128( t2, r6 );
+
+				t5 = _mm_or_si128( t3, t4 );
+				// Store 8 FP16 values, the destination is aligned
+				_mm_store_si128( ( __m128i* )rdi, t5 );
+				rdi += destStride;
+			}
+		}
+	}
+}
+
+// At least for the hybrid decoder, this method absolutely dominates the CPU time.
+// And not due to the integer shuffles - the bottleneck is loading data from the source matrix.
+HRESULT MulMatBase::transposePanelAvx2( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	assert( stridesA[ 0 ] == 1 );
+
+	const size_t heightFloats = (size_t)panelHeightRegisters * 8;
+	i *= heightFloats;
+
+	const uint16_t* rsi = (const uint16_t*)pa;
+	rsi += m3 * stridesA[ 3 ];
+	rsi += m2 * stridesA[ 2 ];
+	rsi += i * stridesA[ 1 ];
+
+	const size_t resultStride = heightFloats;
+
+	if( i + heightFloats <= resultSize[ 0 ] )
+	{
+		// A complete panel
+		for( size_t i = 0; i < panelHeightRegisters; i++ )
+		{
+			transpose8Avx2( rdi, length, rsi, stridesA[ 1 ], resultStride );
+			// Advance by 8 floats in the output buffer
+			rdi += 8;
+			// Advance by 8 rows in the source matrix
+			rsi += 8 * stridesA[ 1 ];
+		}
+	}
+	else
+	{
+		// A partial panel, at the bottom of the first argument matrix
+		const size_t remainder = resultSize[ 0 ] - i;
+		assert( remainder > 0 && remainder < heightFloats );
+		zeroAlignedMemory( rdi, resultStride * length * sizeof( uint16_t ) );
+
+		const size_t completePanels = remainder / 8;
+		for( size_t i = 0; i < completePanels; i++ )
+		{
+			transpose8Avx2( rdi, length, rsi, stridesA[ 1 ], resultStride );
+			rdi += 8;
+			rsi += 8 * stridesA[ 1 ];
+		}
+		const size_t lastPanel = remainder % 8;
+		if( 0 != lastPanel )
+			transpose8PartialAvx2( rdi, length, lastPanel, rsi, stridesA[ 1 ], resultStride );
+	}
+	return S_OK;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/mulMatImpl.cpp b/Whisper/CPU/mulMatImpl.cpp
new file mode 100644
index 0000000..fc50b03
--- /dev/null
+++ b/Whisper/CPU/mulMatImpl.cpp
@@ -0,0 +1,213 @@
+#include "stdafx.h"
+#include <intrin.h>
+#include "mulMatImpl.h"
+#include "mulMat.kernel.hpp"
+
+#define DBG_TRACK_TEMPLATE_INSTANTIATION 0
+
+#if DBG_TRACK_TEMPLATE_INSTANTIATION
+#include <unordered_set>
+static std::unordered_set<uint16_t> g_mulMatTemplates;
+#endif
+
+namespace
+{
+	using namespace CpuCompute;
+
+	bool checkAvx2Support()
+	{
+		int cpuInfo[ 4 ];
+		__cpuid( cpuInfo, 7 );
+		return ( cpuInfo[ 1 ] & ( 1 << 5 ) ) != 0;
+	}
+
+	// a / b, rounded up to the next integer
+	inline uint32_t divRoundUp( uint32_t a, uint32_t b )
+	{
+		assert( b != 0 );
+		return ( a + ( b - 1 ) ) / b;
+	}
+}
+
+const bool MulMatBase::haveAvx2 = checkAvx2Support();
+
+MulMatBase::MulMatBase( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor, uint8_t panelHeightRegs, uint8_t tileWidthFloats ) :
+	resultPointer( result.fp32() ),
+	pa( a.data() ),
+	pb( b.data() ),
+	runner( pfor )
+{
+	length = a.ne[ 0 ];
+	resultStrides[ 0 ] = result.nb[ 1 ];
+	resultStrides[ 1 ] = result.nb[ 2 ];
+	resultStrides[ 2 ] = result.nb[ 3 ];
+	store( resultSize, result.sizeVec() );
+	store( stridesA, a.stridesVec() );
+	store( stridesB, b.stridesVec() );
+
+	countPanels = divRoundUp( resultSize[ 0 ], panelHeightRegs * 8 );
+	completeTilesPerPanel = resultSize[ 1 ] / tileWidthFloats;
+	lastColumnsInPanel = (uint8_t)( resultSize[ 1 ] % tileWidthFloats );
+	this->panelHeightRegisters = panelHeightRegs;
+	this->tileWidth = tileWidthFloats;
+
+	// Pick a method which reshapes a panel of the matrix A into the shape we need to compute the product
+	// Store the pointer to that method in the field of this class
+	if( a.nb[ 0 ] == 1 )
+	{
+		if( haveAvx2 )
+			pfnMakePanel = &MulMatBase::transposePanelAvx2;
+		else
+			pfnMakePanel = &MulMatBase::transposePanel;
+	}
+	else if( a.nb[ 1 ] == 1 )
+	{
+		switch( panelHeightRegs )
+		{
+		case 1:
+			pfnMakePanel = &MulMatBase::copyPanelColumnMajor8;
+			break;
+		case 2:
+			pfnMakePanel = &MulMatBase::copyPanelColumnMajor16;
+			break;
+		case 4:
+			pfnMakePanel = &MulMatBase::copyPanelColumnMajor32;
+			break;
+		default:
+			throw E_NOTIMPL;
+		}
+	}
+	else
+		pfnMakePanel = &MulMatBase::gatherPanel;
+
+	// That last version is generic and very simple, unlikely to have weird bugs
+	// pfnMakePanel = &MulMatBase::gatherPanel;
+
+#if DBG_TRACK_TEMPLATE_INSTANTIATION
+	uint16_t key = panelHeightRegs;
+	key = key << 8;
+	key |= tileWidthFloats;
+	if( !g_mulMatTemplates.emplace( key ).second )
+		return;
+	logDebug( u8"MulMatImpl<panelHeightRegs = %i, tileWidthFloats = %i>", (int)panelHeightRegs, (int)tileWidthFloats );
+#endif
+}
+
+HRESULT MulMatBase::run( ParallelForRunner& pfor )
+{
+	size_t length = (size_t)countPanels * resultSize[ 2 ] * resultSize[ 3 ];
+	return pfor.parallelFor( *this, length );
+}
+
+const float* MulMatBase::getLayerB( size_t m2, size_t m3 ) const
+{
+	const float* rsi = (const float*)this->pb;
+	rsi += m2 * stridesB[ 2 ];
+	rsi += m3 * stridesB[ 3 ];
+	return rsi;
+}
+
+// This method is the main one, it’s called by the thread pool
+template<uint8_t panelHeightRegs, uint8_t tileWidthFloats>
+HRESULT __stdcall MulMatImpl<panelHeightRegs, tileWidthFloats>::compute( size_t i, size_t end ) const noexcept
+{
+	// Allocate a thread-local buffer for the transposed panel
+	constexpr size_t panelHeightFloats = panelHeightRegs * 8;
+	uint16_t* const panel = (uint16_t*)runner.threadLocalBuffer( floatsPerPanel() * 2 );
+	const size_t resultStride = resultStrides[ 0 ];
+
+	// Load a few numbers from this class into local variables, while upcasting from DWORD into size_t
+	const size_t length = this->length;
+	const std::array<size_t, 2> stridesB{ this->stridesB[ 0 ], this->stridesB[ 1 ] };
+
+	// This outer loop iterates over the panels assigned to the current thread
+	// For example, matrix A of size [ 1024, 1024 ] may be split into panels of size [ 1024, 16 ]
+	// Each iteration of that loop computes matrix product of that panel, with the complete matrix B
+	for( ; i < end; i++ )
+	{
+		const size_t iPanel = i % countPanels;
+		size_t j = i / countPanels;
+		const size_t m2 = j % (size_t)resultSize[ 2 ];
+		const size_t m3 = j / (size_t)resultSize[ 2 ];
+
+		CHECK( ( this->*pfnMakePanel )( panel, iPanel, m2, m3 ) );
+		// We got a column-major panel in the thread local buffer, of size [ length, panelHeightRegs * 8 ]
+		// Hopefully, these buffers should all fit at least in L3 cache
+		// The longest matrix I saw in the debugger had 4096 elements, with panelHeightRegs = 4 that's 256 kb of data in the panel
+		const float* pb = getLayerB( m2, m3 );
+		float* rdi = getPanelDest( iPanel, m2, m3 );
+
+		const size_t storeWidth = std::min( panelHeightFloats, (size_t)resultSize[ 0 ] - iPanel * panelHeightFloats );
+		std::array<__m256, panelHeightRegs> vecPanel;
+#if 1
+		ResultTile<panelHeightRegs, tileWidthFloats> tile;
+
+		// This loop iterates over tiles within the panel.
+		// Each iteration of the loop computes an output tile of the result matrix.
+		for( j = 0; j < completeTilesPerPanel; j++, pb += tileWidthFloats * stridesB[ 1 ], rdi += resultStride * tileWidthFloats )
+		{
+			setZero( tile.arr );
+			const uint16_t* rsiA = panel;
+			const uint16_t* const rsiAEnd = panel + length * panelHeightFloats;
+			const float* rsiB = pb;
+			// This loop runs for `length` iterations, iterates over the first dimensions of both matrices, accumulating these dot products we're after
+			for( ; rsiA < rsiAEnd; rsiA += panelHeightFloats, rsiB += stridesB[ 0 ] )
+			{
+				loadPanel( rsiA, vecPanel );
+				tile.kernel( vecPanel, rsiB, stridesB[ 1 ] );
+			}
+			tile.store( rdi, storeWidth, tileWidthFloats, resultStride );
+		}
+
+		if( 0 != lastColumnsInPanel )
+		{
+			setZero( tile.arr );
+			const uint16_t* rsiA = panel;
+			const uint16_t* rsiAEnd = panel + length * panelHeightFloats;
+			const float* rsiB = pb;
+			for( ; rsiA < rsiAEnd; rsiA += panelHeightFloats, rsiB += stridesB[ 0 ] )
+			{
+				loadPanel( rsiA, vecPanel );
+				tile.kernelPartial( vecPanel, rsiB, stridesB[ 1 ], lastColumnsInPanel );
+			}
+			tile.store( rdi, storeWidth, lastColumnsInPanel, resultStride );
+		}
+#else
+		// This version bypasses horizontal tiling, instead implements a brute force algorithm to multiply the current panel by the complete B matrix
+		// Not terribly efficient, only implemented for debugging purposes
+		const size_t resHeight = resultSize[ 1 ];
+		std::array<__m256, panelHeightRegs> tile;
+		for( size_t j = 0; j < resHeight; j++, pb += stridesB[ 1 ], rdi += resultStride )
+		{
+			setZero( tile );
+
+			const uint16_t* rsiA = panel;
+			const uint16_t* const rsiAEnd = panel + length * panelHeightFloats;
+			const float* rsiB = pb;
+			for( size_t k = 0; k < length; k++, rsiA += panelHeightFloats, rsiB += stridesB[ 0 ] )
+			{
+				loadPanel( rsiA, vecPanel );
+				const __m256 b = _mm256_broadcast_ss( rsiB );
+				for( size_t r = 0; r < panelHeightRegs; r++ )
+					tile[ r ] = _mm256_fmadd_ps( vecPanel[ r ], b, tile[ r ] );
+			}
+
+			alignas( 32 ) std::array<float, panelHeightFloats> arr;
+			for( size_t k = 0; k < panelHeightRegs; k++ )
+				_mm256_store_ps( &arr[ k * 8 ], tile[ k ] );
+			memcpy( rdi, arr.data(), storeWidth * 4 );
+		}
+#endif
+	}
+	return S_OK;
+}
+
+// Instantiate the templates we need
+template class MulMatImpl<4, 1>;
+template class MulMatImpl<1, 1>;
+template class MulMatImpl<4, 2>;
+template class MulMatImpl<1, 2>;
+template class MulMatImpl<2, 3>;
+template class MulMatImpl<1, 3>;
+template class MulMatImpl<2, 4>;
+template class MulMatImpl<1, 4>;
\ No newline at end of file
diff --git a/Whisper/CPU/mulMatImpl.h b/Whisper/CPU/mulMatImpl.h
new file mode 100644
index 0000000..8e0062f
--- /dev/null
+++ b/Whisper/CPU/mulMatImpl.h
@@ -0,0 +1,106 @@
+#pragma once
+// Matrix*matrix multiplication is the most expensive algorithm in the model, by far.
+// For this reason, the code in this source file, and in the mulMat.kernel.hpp header, is optimized for performance. Readability suffers.
+// The implementation is inspired by following two articles:
+// https://gist.github.com/nadavrot/5b35d44e8ba3dd718e595e40184d03f0
+// https://link.springer.com/article/10.1007/s11227-022-05003-3
+#include "ParallelForRunner.h"
+#include "Tensor.h"
+
+namespace CpuCompute
+{
+	// Abstract base class for all implementations, to reduce binary size
+	class MulMatBase : public iComputeRange
+	{
+	protected:
+		// Pointers to the payload of the output matrix
+		float* const resultPointer;
+
+		// Lengths of the dot products to compute, equal to width of both source matrices
+		uint32_t length;
+
+		// Last 3 strides of the output matrix, expressed as count of elements. The first one is always 1 because the output matrix is continuous.
+		std::array<uint32_t, 3> resultStrides;
+
+		// Size of the output matrix
+		std::array<uint32_t, 4> resultSize;
+
+		// Pointers to the payload of the source matrices
+		const void* const pa;
+		const void* const pb;
+
+		// Matrix strides, expressed as count of elements
+		std::array<uint32_t, 4> stridesA, stridesB;
+
+		// Total count of panels in the layer of the output matrix.
+		// The last panel might be incomplete, with smaller height.
+		// The thread-local buffer however is always complete, unused elements will be zeros.
+		uint32_t countPanels;
+
+		// Complete tiles in the length of the panel
+		uint32_t completeTilesPerPanel;
+
+		// Count of the last remainder columns in the panel, can be 0
+		uint8_t lastColumnsInPanel;
+
+		// Same as panelHeightRegs template argument - height of the panels, in AVX vectors
+		uint8_t panelHeightRegisters;
+
+		// Same as tileWidthFloats template argument - width of the tile, in floats
+		uint8_t tileWidth;
+
+		// Method pointer to reshape a panel from the source matrix into a thread-local buffer
+		using pfnTransposePanel = HRESULT( MulMatBase::* )( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		pfnTransposePanel pfnMakePanel;
+		// The object which implements multithreading for this job, and supplies memory for thread-local buffers
+		ParallelForRunner& runner;
+
+		// Count of FP16 values in the thread-local panel buffer
+		uint32_t floatsPerPanel() const
+		{
+			return length * panelHeightRegisters * 8;
+		}
+
+		// Transpose a horizontal panel of the first matrix, when the rows are continuous in that matrix
+		HRESULT transposePanel( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		HRESULT transposePanelAvx2( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		// Copy a horizontal panel of the first matrix without transpose, for column major layout of that matrix
+		HRESULT copyPanelColumnMajor8( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		HRESULT copyPanelColumnMajor16( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		HRESULT copyPanelColumnMajor32( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+		// Transpose a panel of the first matrix for irregular layout of that matrix, when neither rows nor columns are at sequential addresses.
+		// This one ain't implemented yet.
+		HRESULT gatherPanel( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const;
+
+		const uint16_t* getPanelA( size_t i, size_t m2, size_t m3 ) const;
+		// Pointer to the first element of the second source matrix in the specified layer
+		const float* getLayerB( size_t m2, size_t m3 ) const;
+
+		// Pointer to the first element of the output tile of the result matrix
+		float* getPanelDest( size_t i, size_t m2, size_t m3 ) const
+		{
+			float* rdi = resultPointer;
+			rdi += m2 * resultStrides[ 1 ];
+			rdi += m3 * resultStrides[ 2 ];
+			rdi += i * panelHeightRegisters * 8;
+			return rdi;
+		}
+
+		static const bool haveAvx2;
+	public:
+		MulMatBase( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor, uint8_t panelHeightRegs, uint8_t tileWidthFloats );
+		HRESULT run( ParallelForRunner& pfor );
+	};
+
+	// This class actually contains the kernels implementations
+	template<uint8_t panelHeightRegs, uint8_t tileWidthFloats>
+	class MulMatImpl : public MulMatBase
+	{
+		HRESULT __stdcall compute( size_t i, size_t end ) const noexcept override final;
+
+	public:
+		MulMatImpl( Tensor& result, const Tensor& a, const Tensor& b, ParallelForRunner& pfor ) :
+			MulMatBase( result, a, b, pfor, panelHeightRegs, tileWidthFloats )
+		{ }
+	};
+}
\ No newline at end of file
diff --git a/Whisper/CPU/mulMatImpl.panel.cpp b/Whisper/CPU/mulMatImpl.panel.cpp
new file mode 100644
index 0000000..f3baf21
--- /dev/null
+++ b/Whisper/CPU/mulMatImpl.panel.cpp
@@ -0,0 +1,274 @@
+#include "stdafx.h"
+#include <intrin.h>
+#include "mulMatImpl.h"
+#include "mulMatUtils.hpp"
+using namespace CpuCompute;
+
+// We want to keep code size reasonable, that's why these panel reshaping methods are in the base class
+HRESULT MulMatBase::transposePanel( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	assert( stridesA[ 0 ] == 1 );
+
+	const size_t heightFloats = (size_t)panelHeightRegisters * 8;
+	i *= heightFloats;
+
+	const uint16_t* rsi = (const uint16_t*)pa;
+	rsi += m3 * stridesA[ 3 ];
+	rsi += m2 * stridesA[ 2 ];
+	rsi += i * stridesA[ 1 ];
+
+	const size_t resultStride = heightFloats;
+
+	if( i + heightFloats <= resultSize[ 0 ] )
+	{
+		// A complete panel
+		for( size_t i = 0; i < panelHeightRegisters; i++ )
+		{
+			transpose8( rdi, length, rsi, stridesA[ 1 ], resultStride );
+			// Advance by 8 floats in the output buffer
+			rdi += 8;
+			// Advance by 8 rows in the source matrix
+			rsi += 8 * stridesA[ 1 ];
+		}
+	}
+	else
+	{
+		// A partial panel, at the bottom of the first argument matrix
+		const size_t remainder = resultSize[ 0 ] - i;
+		assert( remainder > 0 && remainder < heightFloats );
+		zeroAlignedMemory( rdi, resultStride * length * sizeof( uint16_t ) );
+
+		const size_t completePanels = remainder / 8;
+		for( size_t i = 0; i < completePanels; i++ )
+		{
+			transpose8( rdi, length, rsi, stridesA[ 1 ], resultStride );
+			rdi += 8;
+			rsi += 8 * stridesA[ 1 ];
+		}
+		const size_t lastPanel = remainder % 8;
+		if( 0 != lastPanel )
+			transpose8Partial( rdi, length, lastPanel, rsi, stridesA[ 1 ], resultStride );
+	}
+	return S_OK;
+}
+
+inline const uint16_t* MulMatBase::getPanelA( size_t i, size_t m2, size_t m3 ) const
+{
+	const uint16_t* rsi = (const uint16_t*)pa;
+	rsi += m3 * stridesA[ 3 ];
+	rsi += m2 * stridesA[ 2 ];
+	rsi += i * stridesA[ 1 ];
+	return rsi;
+}
+
+HRESULT MulMatBase::copyPanelColumnMajor8( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	assert( stridesA[ 1 ] == 1 );
+	assert( panelHeightRegisters == 1 );
+
+	constexpr size_t heightFloats = 8;
+	i *= heightFloats;
+	const uint16_t* rsi = getPanelA( i, m2, m3 );
+
+	constexpr size_t resultStride = heightFloats;
+
+	if( i + heightFloats <= resultSize[ 0 ] )
+	{
+		// A complete panel, height = 8 elements
+		copyColumnMajor( rdi, length, rsi, stridesA[ 0 ], resultStride );
+	}
+	else
+	{
+		// A partial panel, at the bottom of the first argument matrix
+		const size_t remainder = resultSize[ 0 ] - i;
+		assert( remainder > 0 && remainder < heightFloats );
+		copyColumnMajorPartial( rdi, length, remainder, rsi, stridesA[ 0 ], resultStride );
+	}
+	return S_OK;
+}
+
+__forceinline __m128i load8Partial( const uint16_t* x, size_t len )
+{
+	assert( len > 0 && len < 8 );
+	__m128i ix = _mm_setzero_si128();
+	switch( len )
+	{
+	case 1: // load 2 bytes
+		ix = _mm_cvtsi32_si128( *x );
+		break;
+	case 2: // load 4 bytes
+		ix = _mm_cvtsi32_si128( *(const int*)x );
+		break;
+	case 3: // load 6 bytes
+		ix = _mm_cvtsi32_si128( *(const int*)x );
+		ix = _mm_insert_epi16( ix, x[ 2 ], 2 );
+		break;
+	case 4: // load 8 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		break;
+	case 5: // load 10 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi16( ix, x[ 4 ], 4 );
+		break;
+	case 6: // load 12 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+		break;
+	case 7: // load 14 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+		ix = _mm_insert_epi16( ix, x[ 6 ], 6 );
+		break;
+	}
+	return ix;
+}
+
+__forceinline __m256i load16Partial( const uint16_t* rsi, size_t len )
+{
+	assert( len > 0 && len < 16 );
+
+	if( len < 8 )
+	{
+		__m128i low = load8Partial( rsi, len );
+		return _mm256_setr_m128i( low, _mm_setzero_si128() );
+	}
+	else if( len > 8 )
+	{
+		__m128i low = load16( (const int*)rsi );
+		__m128i high = load8Partial( rsi + 8, len - 8 );
+		return _mm256_setr_m128i( low, high );
+	}
+	else
+	{
+		__m128i low = load16( (const int*)rsi );
+		return _mm256_setr_m128i( low, _mm_setzero_si128() );
+	}
+}
+
+HRESULT MulMatBase::copyPanelColumnMajor16( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	assert( stridesA[ 1 ] == 1 );
+	assert( panelHeightRegisters == 2 );
+
+	constexpr size_t heightFloats = 16;
+	i *= heightFloats;
+
+	const uint16_t* rsi = getPanelA( i, m2, m3 );
+	uint16_t* const rdiEnd = rdi + 16 * length;
+
+	if( i + heightFloats <= resultSize[ 0 ] )
+	{
+		// A complete panel, height = 16 elements
+		for( ; rdi < rdiEnd; rdi += 16, rsi += stridesA[ 0 ] )
+		{
+			__m256i v = _mm256_loadu_si256( ( const __m256i* )rsi );
+			_mm256_store_si256( ( __m256i* )rdi, v );
+		}
+	}
+	else
+	{
+		// A partial panel, at the bottom of the first argument matrix
+		const size_t remainder = resultSize[ 0 ] - i;
+		assert( remainder > 0 && remainder < heightFloats );
+
+		for( ; rdi < rdiEnd; rdi += 16, rsi += stridesA[ 0 ] )
+		{
+			__m256i v = load16Partial( rsi, remainder );
+			_mm256_store_si256( ( __m256i* )rdi, v );
+		}
+	}
+	return S_OK;
+}
+
+HRESULT MulMatBase::copyPanelColumnMajor32( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	assert( stridesA[ 1 ] == 1 );
+	assert( panelHeightRegisters == 4 );
+
+	constexpr size_t heightFloats = 32;
+	i *= heightFloats;
+
+	const uint16_t* rsi = getPanelA( i, m2, m3 );
+	uint16_t* const rdiEnd = rdi + 32 * length;
+
+	if( i + heightFloats <= resultSize[ 0 ] )
+	{
+		// A complete panel, height = 32 elements
+		for( ; rdi < rdiEnd; rdi += 32, rsi += stridesA[ 0 ] )
+		{
+			__m256i v = _mm256_loadu_si256( ( const __m256i* )rsi );
+			_mm256_store_si256( ( __m256i* )rdi, v );
+			v = _mm256_loadu_si256( ( const __m256i* )( rsi + 16 ) );
+			_mm256_store_si256( ( __m256i* )( rdi + 16 ), v );
+		}
+	}
+	else
+	{
+		// A partial panel, at the bottom of the first argument matrix
+		const size_t remainder = resultSize[ 0 ] - i;
+		assert( remainder > 0 && remainder < heightFloats );
+
+		// _mm256_setzero_si256 probably compiles into vpxor, that's AVX2, we don't want that here
+		const __m256 zero = _mm256_setzero_ps();
+
+		for( ; rdi < rdiEnd; rdi += 32, rsi += stridesA[ 0 ] )
+		{
+			if( remainder < 16 )
+			{
+				__m256i v = load16Partial( rsi, remainder );
+				_mm256_store_si256( ( __m256i* )rdi, v );
+				_mm256_store_ps( (float*)( rdi + 16 ), zero );
+			}
+			else if( remainder > 16 )
+			{
+				__m256i v = _mm256_loadu_si256( ( const __m256i* )rsi );
+				_mm256_store_si256( ( __m256i* )rdi, v );
+				v = load16Partial( rsi + 16, remainder - 16 );
+				_mm256_store_si256( ( __m256i* )( rdi + 16 ), v );
+			}
+			else
+			{
+				__m256i v = _mm256_loadu_si256( ( const __m256i* )rsi );
+				_mm256_store_si256( ( __m256i* )rdi, v );
+				_mm256_store_ps( (float*)( rdi + 16 ), zero );
+			}
+		}
+	}
+	return S_OK;
+}
+
+HRESULT MulMatBase::gatherPanel( uint16_t* rdi, size_t i, size_t m2, size_t m3 ) const
+{
+	// BTW, I never saw this method called.
+	const size_t heightFloats = (size_t)panelHeightRegisters * 8;
+	const size_t length = this->length;
+
+	zeroAlignedMemory( rdi, length * heightFloats * sizeof( uint16_t ) );
+
+	const size_t height = std::min( heightFloats, resultSize[ 0 ] - i );
+	const size_t strideElement = stridesA[ 0 ];
+	const size_t strideRow = stridesA[ 1 ];
+	const uint16_t* rsi = getPanelA( i * heightFloats, m2, m3 );
+
+	if( strideElement < strideRow )
+	{
+		for( size_t r = 0; r < height; r++, rsi += strideRow, rdi++ )
+		{
+			const uint16_t* sourceRow = rsi;
+			uint16_t* destRow = rdi;
+			for( size_t c = 0; c < length; c++, sourceRow += strideElement, destRow += heightFloats )
+				*destRow = *sourceRow;
+		}
+	}
+	else
+	{
+		for( size_t c = 0; c < length; c++, rsi += strideElement, rdi += heightFloats )
+		{
+			const uint16_t* sourceCol = rsi;
+			uint16_t* destCol = rdi;
+			for( size_t r = 0; r < height; r++, sourceCol += strideRow, destCol++ )
+				*destCol = *sourceCol;
+		}
+	}
+	return S_OK;
+}
\ No newline at end of file
diff --git a/Whisper/CPU/mulMatUtils.hpp b/Whisper/CPU/mulMatUtils.hpp
new file mode 100644
index 0000000..1276323
--- /dev/null
+++ b/Whisper/CPU/mulMatUtils.hpp
@@ -0,0 +1,301 @@
+#pragma once
+#include <immintrin.h>
+#include <stdint.h>
+#include <assert.h>
+
+__forceinline __m128i f16Load( const uint16_t* rsi )
+{
+	return _mm_loadu_si128( ( const __m128i* )rsi );
+}
+
+constexpr size_t maskAlign8 = ~(size_t)7;
+
+__forceinline void transpose8( uint16_t* rdi, size_t w, const uint16_t* rsi, size_t sourceStride, size_t destStride )
+{
+	assert( 0 == ( (size_t)rdi ) % 16 );
+	assert( 0 == destStride % 8 );
+	assert( w <= sourceStride );
+
+	const uint16_t* const rsiEndAligned = rsi + ( w & maskAlign8 );
+	const uint16_t* rsi5 = rsi + sourceStride * 5;
+	uint16_t* rdi5 = rdi + destStride * 5;
+	const size_t rem = w % 8;
+	for( ; rsi < rsiEndAligned; rsi += 8, rsi5 += 8, rdi += 8 * destStride, rdi5 += 8 * destStride )
+	{
+		// Load 8x8 block into 8 registers
+		__m128i r0 = f16Load( rsi );                     // 00, 01, 02, 03, 04, 05, 06, 07
+		__m128i r1 = f16Load( rsi + sourceStride );      // 10, 11, 12, 13, 14, 15, 16, 17
+		__m128i r2 = f16Load( rsi + sourceStride * 2 );  // 20, 21, 22, 23, 24, 25, 26, 27
+		__m128i r3 = f16Load( rsi5 - sourceStride * 2 ); // 30, 31, 32, 33, 34, 35, 36, 37
+		__m128i r4 = f16Load( rsi5 - sourceStride );     // 40, 41, 42, 43, 44, 45, 46, 47
+		__m128i r5 = f16Load( rsi5 );                    // 50, 51, 52, 53, 54, 55, 56, 57
+		__m128i r6 = f16Load( rsi5 + sourceStride );     // 60, 61, 62, 63, 64, 65, 66, 67
+		__m128i r7 = f16Load( rsi5 + sourceStride * 2 ); // 70, 71, 72, 73, 74, 75, 76, 77
+
+		// Transpose FP16 values in registers
+		__m128i t0 = _mm_unpacklo_epi16( r0, r1 ); // 00, 10, 01, 11, 02, 12, 03, 13
+		__m128i t1 = _mm_unpackhi_epi16( r0, r1 ); // 04, 14, 05, 15, 06, 16, 07, 17
+		__m128i t2 = _mm_unpacklo_epi16( r2, r3 ); // 20, 30, 21, 31, 22, 32, 23, 33
+		__m128i t3 = _mm_unpackhi_epi16( r2, r3 ); // 24, 34, 25, 35, 26, 36, 27, 37
+		__m128i t4 = _mm_unpacklo_epi16( r4, r5 ); // 40, 50, 41, 52, 42, 52, 43, 53
+		__m128i t5 = _mm_unpackhi_epi16( r4, r5 ); // 44, 54, 45, 55, 46, 56, 47, 57
+		__m128i t6 = _mm_unpacklo_epi16( r6, r7 ); // 60, 70, 61, 71, 62, 72, 63, 73
+		__m128i t7 = _mm_unpackhi_epi16( r6, r7 ); // 64, 74, 65, 75, 66, 76, 67, 77
+
+		r0 = _mm_unpacklo_epi32( t0, t2 ); // 00, 10, 20, 30, 01, 11, 21, 31
+		r1 = _mm_unpackhi_epi32( t0, t2 ); // 02, 12, 22, 32, 03, 13, 23, 33
+		r2 = _mm_unpacklo_epi32( t1, t3 ); // 04, 14, 24, 34, 05, 15, 25, 35
+		r3 = _mm_unpackhi_epi32( t1, t3 ); // 06, 16, 26, 36, 07, 17, 27, 37
+		r4 = _mm_unpacklo_epi32( t4, t6 ); // 40, 50, 60, 70, 41, 51, 61, 71
+		r5 = _mm_unpackhi_epi32( t4, t6 ); // 42, 52, 62, 72, 43, 53, 63, 73
+		r6 = _mm_unpacklo_epi32( t5, t7 ); // 44, 54, 64, 74, 45, 55, 65, 75
+		r7 = _mm_unpackhi_epi32( t5, t7 ); // 46, 56, 66, 76, 47, 57, 67, 77
+
+		t0 = _mm_unpacklo_epi64( r0, r4 ); // 00, 10, 20, 30, 40, 50, 60, 70
+		t1 = _mm_unpackhi_epi64( r0, r4 ); // 01, 11, 21, 31, 41, 52, 61, 71
+		t2 = _mm_unpacklo_epi64( r1, r5 ); // 02, 12, 22, 32, 42, 52, 62, 72
+		t3 = _mm_unpackhi_epi64( r1, r5 ); // 03, 13, 23, 33, 43, 53, 63, 73
+		t4 = _mm_unpacklo_epi64( r2, r6 );
+		t5 = _mm_unpackhi_epi64( r2, r6 );
+		t6 = _mm_unpacklo_epi64( r3, r7 );
+		t7 = _mm_unpackhi_epi64( r3, r7 );
+
+		// Store
+		store16( rdi, t0 );
+		store16( rdi + destStride, t1 );
+		store16( rdi + destStride * 2, t2 );
+		store16( rdi5 - destStride * 2, t3 );
+		store16( rdi5 - destStride, t4 );
+		store16( rdi5, t5 );
+		store16( rdi5 + destStride, t6 );
+		store16( rdi5 + destStride * 2, t7 );
+	}
+
+#pragma loop( no_vector )
+	for( size_t i = 0; i < rem; rsi++, rsi5++, rdi += destStride )
+	{
+		const int16_t* p0 = (const int16_t*)rsi;
+		const int16_t* p5 = (const int16_t*)rsi5;
+		// Load a complete column into a vector
+		__m128i v = _mm_cvtsi32_si128( *rsi );
+		v = _mm_insert_epi16( v, *( p0 + sourceStride ), 1 );
+		v = _mm_insert_epi16( v, *( p0 + sourceStride * 2 ), 2 );
+		v = _mm_insert_epi16( v, *( p5 - sourceStride * 2 ), 3 );
+		v = _mm_insert_epi16( v, *( p5 - sourceStride ), 4 );
+		v = _mm_insert_epi16( v, *( p5 ), 5 );
+		v = _mm_insert_epi16( v, *( p5 + sourceStride ), 6 );
+		v = _mm_insert_epi16( v, *( p5 + sourceStride * 2 ), 7 );
+		// Store 8 FP16 values
+		store16( rdi, v );
+	}
+}
+
+inline void transpose8Partial( uint16_t* rdi, size_t w, size_t h, const uint16_t* rsi, size_t sourceStride, size_t destStride )
+{
+	assert( 0 == ( (size_t)rdi ) % 16 );
+	assert( 0 == destStride % 8 );
+	assert( w <= sourceStride );
+	assert( h > 0 && h < 8 );
+
+	const uint16_t* const rsiEndAligned = rsi + ( w & maskAlign8 );
+	const uint16_t* rsi5 = rsi + sourceStride * 5;
+	uint16_t* rdi5 = rdi + destStride * 5;
+	const size_t rem = w % 8;
+	for( ; rsi < rsiEndAligned; rsi += 8, rsi5 += 8, rdi += 8 * destStride, rdi5 += 8 * destStride )
+	{
+		// Load the block into 8 registers, set unused rows to zero
+		__m128i r0 = f16Load( rsi );
+		__m128i r1 = _mm_setzero_si128();
+		__m128i r2 = _mm_setzero_si128();
+		__m128i r3 = _mm_setzero_si128();
+		__m128i r4 = _mm_setzero_si128();
+		__m128i r5 = _mm_setzero_si128();
+		__m128i r6 = _mm_setzero_si128();
+		// These branches, whether direct or indirect, are very predictable: same outcome for all iterations of the outer loop
+		switch( h )
+		{
+		case 7:
+			r6 = f16Load( rsi5 + sourceStride );
+		case 6:
+			r5 = f16Load( rsi5 );
+		case 5:
+			r4 = f16Load( rsi5 - sourceStride );
+		case 4:
+			r3 = f16Load( rsi5 - sourceStride * 2 );
+		case 3:
+			r2 = f16Load( rsi + sourceStride * 2 );
+		case 2:
+			r1 = f16Load( rsi + sourceStride );
+		}
+		__m128i r7 = _mm_setzero_si128();
+
+		// Transpose FP16 values in registers
+		__m128i t0 = _mm_unpacklo_epi16( r0, r1 ); // 00, 10, 01, 11, 02, 12, 03, 13
+		__m128i t1 = _mm_unpackhi_epi16( r0, r1 ); // 04, 14, 05, 15, 06, 16, 07, 17
+		__m128i t2 = _mm_unpacklo_epi16( r2, r3 ); // 20, 30, 21, 31, 22, 32, 23, 33
+		__m128i t3 = _mm_unpackhi_epi16( r2, r3 ); // 24, 34, 25, 35, 26, 36, 27, 37
+		__m128i t4 = _mm_unpacklo_epi16( r4, r5 ); // 40, 50, 41, 52, 42, 52, 43, 53
+		__m128i t5 = _mm_unpackhi_epi16( r4, r5 ); // 44, 54, 45, 55, 46, 56, 47, 57
+		__m128i t6 = _mm_unpacklo_epi16( r6, r7 ); // 60, 70, 61, 71, 62, 72, 63, 73
+		__m128i t7 = _mm_unpackhi_epi16( r6, r7 ); // 64, 74, 65, 75, 66, 76, 67, 77
+
+		r0 = _mm_unpacklo_epi32( t0, t2 ); // 00, 10, 20, 30, 01, 11, 21, 31
+		r1 = _mm_unpackhi_epi32( t0, t2 ); // 02, 12, 22, 32, 03, 13, 23, 33
+		r2 = _mm_unpacklo_epi32( t1, t3 ); // 04, 14, 24, 34, 05, 15, 25, 35
+		r3 = _mm_unpackhi_epi32( t1, t3 ); // 06, 16, 26, 36, 07, 17, 27, 37
+		r4 = _mm_unpacklo_epi32( t4, t6 ); // 40, 50, 60, 70, 41, 51, 61, 71
+		r5 = _mm_unpackhi_epi32( t4, t6 ); // 42, 52, 62, 72, 43, 53, 63, 73
+		r6 = _mm_unpacklo_epi32( t5, t7 ); // 44, 54, 64, 74, 45, 55, 65, 75
+		r7 = _mm_unpackhi_epi32( t5, t7 ); // 46, 56, 66, 76, 47, 57, 67, 77
+
+		t0 = _mm_unpacklo_epi64( r0, r4 ); // 00, 10, 20, 30, 40, 50, 60, 70
+		t1 = _mm_unpackhi_epi64( r0, r4 ); // 01, 11, 21, 31, 41, 52, 61, 71
+		t2 = _mm_unpacklo_epi64( r1, r5 ); // 02, 12, 22, 32, 42, 52, 62, 72
+		t3 = _mm_unpackhi_epi64( r1, r5 ); // 03, 13, 23, 33, 43, 53, 63, 73
+		t4 = _mm_unpacklo_epi64( r2, r6 );
+		t5 = _mm_unpackhi_epi64( r2, r6 );
+		t6 = _mm_unpacklo_epi64( r3, r7 );
+		t7 = _mm_unpackhi_epi64( r3, r7 );
+
+		// Store
+		store16( rdi, t0 );
+		store16( rdi + destStride, t1 );
+		store16( rdi + destStride * 2, t2 );
+		store16( rdi5 - destStride * 2, t3 );
+		store16( rdi5 - destStride, t4 );
+		store16( rdi5, t5 );
+		store16( rdi5 + destStride, t6 );
+		store16( rdi5 + destStride * 2, t7 );
+	}
+
+#pragma loop( no_vector )
+	for( size_t i = 0; i < rem; rsi++, rsi5++, rdi += destStride )
+	{
+		const int16_t* p0 = (const int16_t*)rsi;
+		const int16_t* p5 = (const int16_t*)rsi5;
+		// Load a partial column into vector
+		__m128i v = _mm_cvtsi32_si128( *rsi );
+		switch( h )
+		{
+		case 7:
+			v = _mm_insert_epi16( v, *( p5 + sourceStride ), 6 );
+		case 6:
+			v = _mm_insert_epi16( v, *( p5 ), 5 );
+		case 5:
+			v = _mm_insert_epi16( v, *( p5 - sourceStride ), 4 );
+		case 4:
+			v = _mm_insert_epi16( v, *( p5 - sourceStride * 2 ), 3 );
+		case 3:
+			v = _mm_insert_epi16( v, *( p0 + sourceStride * 2 ), 2 );
+		case 2:
+			v = _mm_insert_epi16( v, *( p0 + sourceStride ), 1 );
+		}
+		// Store 8 FP16 values
+		store16( rdi, v );
+	}
+}
+
+// Same as above, but skip the transpose. The source stride is distance between columns of the matrix.
+__forceinline void copyColumnMajor( uint16_t* rdi, size_t w, const uint16_t* rsi, size_t sourceStride, size_t destStride )
+{
+	assert( 0 == ( (size_t)rdi ) % 16 );
+	assert( 0 == destStride % 8 );
+
+	constexpr size_t maskAlign4 = ~(size_t)3;
+
+	const uint16_t* const rsiEndAligned = rsi + sourceStride * ( w & maskAlign4 );
+	const uint16_t* const rsiEnd = rsi + sourceStride * w;
+	for( ; rsi < rsiEndAligned; rsi += sourceStride * 4, rdi += destStride * 4 )
+	{
+		__m128i c = f16Load( rsi );
+		store16( rdi, c );
+
+		c = f16Load( rsi + sourceStride );
+		store16( rdi + destStride, c );
+
+		c = f16Load( rsi + sourceStride * 2 );
+		store16( rdi + destStride * 2, c );
+
+		c = f16Load( rsi + sourceStride * 3 );
+		store16( rdi + destStride * 3, c );
+	}
+
+	for( ; rsi < rsiEnd; rsi += sourceStride, rdi += destStride )
+	{
+		__m128i c = f16Load( rsi );
+		store16( rdi, c );
+	}
+}
+
+__forceinline __m128i loadPartial( const uint16_t* x, size_t count )
+{
+	assert( count < 8 );
+	__m128i ix;
+	switch( count )
+	{
+	case 1: // load 2 bytes
+		ix = _mm_cvtsi32_si128( *x );
+		break;
+	case 2: // load 4 bytes
+		ix = _mm_cvtsi32_si128( *(const int*)x );
+		break;
+	case 3: // load 6 bytes
+		ix = _mm_cvtsi32_si128( *(const int*)x );
+		ix = _mm_insert_epi16( ix, x[ 2 ], 2 );
+		break;
+	case 4: // load 8 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		break;
+	case 5: // load 10 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi16( ix, x[ 4 ], 4 );
+		break;
+	case 6: // load 12 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+		break;
+	case 7: // load 14 bytes
+		ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+		ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+		ix = _mm_insert_epi16( ix, x[ 6 ], 6 );
+		break;
+	default:
+		return _mm_setzero_si128();
+	}
+	return ix;
+}
+
+inline void copyColumnMajorPartial( uint16_t* rdi, size_t w, size_t h, const uint16_t* rsi, size_t sourceStride, size_t destStride )
+{
+	assert( 0 == ( (size_t)rdi ) % 32 );
+	assert( 0 == destStride % 8 );
+	assert( h > 0 && h < 8 );
+
+	const uint16_t* const rsiEnd = rsi + sourceStride * w;
+	for( ; rsi < rsiEnd; rsi += sourceStride, rdi += destStride )
+	{
+		// Can't use mask loads because loading 2-byte elements
+		// Still, that switch() in loadPartial makes a very predictable branch, same outcome for all iterations of this loop.
+		__m128i c = loadPartial( rsi, h );
+		store16( rdi, c );
+	}
+}
+
+// Store zeros into block of memory, with aligned AVX store instructions
+__forceinline void zeroAlignedMemory( void* pv, size_t cb )
+{
+	assert( 0 == cb % 16 );
+	assert( 0 == ( (size_t)pv % 32 ) );
+
+	uint8_t* rdi = (uint8_t*)pv;
+	constexpr size_t maskAlign32 = ~(size_t)31;
+	uint8_t* const rdiEndAligned = rdi + ( cb & maskAlign32 );
+	uint8_t* const rdiEnd = rdi + cb;
+
+	const __m256 zero = _mm256_setzero_ps();
+	for( ; rdi < rdiEndAligned; rdi += 32 )
+		_mm256_store_ps( (float*)rdi, zero );
+
+	if( rdi < rdiEnd )
+		_mm_store_ps( (float*)rdi, _mm_setzero_ps() );
+}
\ No newline at end of file
diff --git a/Whisper/CPU/simdUtils.cpp b/Whisper/CPU/simdUtils.cpp
new file mode 100644
index 0000000..0e5f77d
--- /dev/null
+++ b/Whisper/CPU/simdUtils.cpp
@@ -0,0 +1,738 @@
+#include "stdafx.h"
+#include "simdUtils.h"
+#include "../ML/LookupTablesData.h"
+#include <cmath>
+#include <memory>
+
+namespace
+{
+	constexpr size_t maskAlign8 = ~(size_t)7;
+
+	__forceinline __m256 load8( const uint16_t* rsi )
+	{
+		__m128i i = _mm_loadu_si128( ( const __m128i* )rsi );
+		return _mm256_cvtph_ps( i );
+	}
+
+	__forceinline void loadPartial( const uint16_t* x, const uint16_t* y, size_t count, __m256& fx, __m256& fy )
+	{
+		assert( count < 8 );
+
+		__m128i ix, iy;
+		switch( count )
+		{
+		case 1: // load 2 bytes
+			ix = _mm_cvtsi32_si128( *x );
+			iy = _mm_cvtsi32_si128( *y );
+			break;
+		case 2: // load 4 bytes
+			ix = _mm_cvtsi32_si128( *(const int*)x );
+			iy = _mm_cvtsi32_si128( *(const int*)y );
+			break;
+		case 3: // load 6 bytes
+			ix = _mm_cvtsi32_si128( *(const int*)x );
+			iy = _mm_cvtsi32_si128( *(const int*)y );
+			ix = _mm_insert_epi16( ix, x[ 2 ], 2 );
+			iy = _mm_insert_epi16( iy, y[ 2 ], 2 );
+			break;
+		case 4: // load 8 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			iy = _mm_cvtsi64_si128( *(const int64_t*)y );
+			break;
+		case 5: // load 10 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			iy = _mm_cvtsi64_si128( *(const int64_t*)y );
+			ix = _mm_insert_epi16( ix, x[ 4 ], 4 );
+			iy = _mm_insert_epi16( iy, y[ 4 ], 4 );
+			break;
+		case 6: // load 12 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			iy = _mm_cvtsi64_si128( *(const int64_t*)y );
+			ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+			iy = _mm_insert_epi32( iy, *(const int*)( y + 4 ), 2 );
+			break;
+		case 7: // load 14 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			iy = _mm_cvtsi64_si128( *(const int64_t*)y );
+			ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+			iy = _mm_insert_epi32( iy, *(const int*)( y + 4 ), 2 );
+			ix = _mm_insert_epi16( ix, x[ 6 ], 6 );
+			iy = _mm_insert_epi16( iy, y[ 6 ], 6 );
+			break;
+		default:
+			fx = fy = _mm256_setzero_ps();
+			return;
+		}
+
+		fx = _mm256_cvtph_ps( ix );
+		fy = _mm256_cvtph_ps( iy );
+	}
+
+	__forceinline __m256 loadPartial( const uint16_t* x, size_t count )
+	{
+		assert( count < 8 );
+		__m128i ix;
+		switch( count )
+		{
+		case 1: // load 2 bytes
+			ix = _mm_cvtsi32_si128( *x );
+			break;
+		case 2: // load 4 bytes
+			ix = _mm_cvtsi32_si128( *(const int*)x );
+			break;
+		case 3: // load 6 bytes
+			ix = _mm_cvtsi32_si128( *(const int*)x );
+			ix = _mm_insert_epi16( ix, x[ 2 ], 2 );
+			break;
+		case 4: // load 8 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			break;
+		case 5: // load 10 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			ix = _mm_insert_epi16( ix, x[ 4 ], 4 );
+			break;
+		case 6: // load 12 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+			break;
+		case 7: // load 14 bytes
+			ix = _mm_cvtsi64_si128( *(const int64_t*)x );
+			ix = _mm_insert_epi32( ix, *(const int*)( x + 4 ), 2 );
+			ix = _mm_insert_epi16( ix, x[ 6 ], 6 );
+			break;
+		default:
+			return _mm256_setzero_ps();
+		}
+		return  _mm256_cvtph_ps( ix );
+	}
+
+	__forceinline __m128 loadFloat2( const float* rsi )
+	{
+		return _mm_castpd_ps( _mm_load_sd( (const double*)rsi ) );
+	}
+	__forceinline __m128 loadFloat3( const float* rsi )
+	{
+		__m128 f = loadFloat2( rsi );
+		f = _mm_insert_ps( f, _mm_load_ss( rsi + 2 ), 0x20 );
+		return f;
+	}
+
+	__forceinline __m256 loadPartial( const float* rsi, size_t count )
+	{
+		assert( count < 8 );
+		__m128 low = _mm_setzero_ps();
+		__m128 high = _mm_setzero_ps();
+		switch( count )
+		{
+		case 1:
+			low = _mm_load_ss( rsi );
+			break;
+		case 2:
+			low = loadFloat2( rsi );
+			break;
+		case 3:
+			low = loadFloat3( rsi );
+			break;
+		case 4:
+			low = _mm_loadu_ps( rsi );
+			break;
+		case 5:
+			low = _mm_loadu_ps( rsi );
+			high = _mm_load_ss( rsi + 4 );
+			break;
+		case 6:
+			low = _mm_loadu_ps( rsi );
+			high = loadFloat2( rsi + 4 );
+			break;
+		case 7:
+			low = _mm_loadu_ps( rsi );
+			high = loadFloat3( rsi + 4 );
+			break;
+		}
+		return _mm256_setr_m128( low, high );
+	}
+
+	__forceinline void storeFloat2( float* rdi, __m128 vec )
+	{
+		_mm_store_sd( (double*)rdi, _mm_castps_pd( vec ) );
+	}
+
+	__forceinline void storePartial( float* rdi, __m256 vec, size_t count )
+	{
+		assert( count < 8 );
+
+		__m128 tmp = _mm256_castps256_ps128( vec );
+		if( count >= 4 )
+		{
+			_mm_storeu_ps( rdi, tmp );
+			if( count == 4 )
+				return;
+			count -= 4;
+			rdi += 4;
+			tmp = _mm256_extractf128_ps( vec, 1 );
+		}
+
+		switch( count )
+		{
+		case 1:
+			_mm_store_ss( rdi, tmp );
+			return;
+		case 2:
+			storeFloat2( rdi, tmp );
+			return;
+		case 3:
+			storeFloat2( rdi, tmp );
+			( (int*)rdi )[ 2 ] = _mm_extract_ps( tmp, 2 );
+			return;
+		}
+	}
+}
+
+void addF16to32( float* rdi, const uint16_t* a, const uint16_t* b, size_t length )
+{
+	const uint16_t* const endAligned = a + ( length & maskAlign8 );
+	const size_t rem = length % 8;
+
+	for( ; a < endAligned; a += 8, b += 8, rdi += 8 )
+	{
+		__m256 f1 = load8( a );
+		__m256 f2 = load8( b );
+		__m256 res = _mm256_add_ps( f1, f2 );
+		_mm256_storeu_ps( rdi, res );
+	}
+
+	if( rem != 0 )
+	{
+		__m256 f1, f2;
+		loadPartial( a, b, rem, f1, f2 );
+		__m256 res = _mm256_add_ps( f1, f2 );
+		storePartial( rdi, res, rem );
+	}
+}
+
+void addF16to32( float* rdi, const uint16_t* a, const float* b, size_t length )
+{
+	const uint16_t* const endAligned = a + ( length & maskAlign8 );
+	const size_t rem = length % 8;
+
+	for( ; a < endAligned; a += 8, b += 8, rdi += 8 )
+	{
+		__m256 f1 = load8( a );
+		__m256 f2 = _mm256_loadu_ps( b );
+		__m256 res = _mm256_add_ps( f1, f2 );
+		_mm256_storeu_ps( rdi, res );
+	}
+
+	if( rem != 0 )
+	{
+		__m256 f1 = loadPartial( a, rem );
+		__m256 f2 = loadPartial( b, rem );
+		__m256 res = _mm256_add_ps( f1, f2 );
+		storePartial( rdi, res, rem );
+	}
+}
+
+alignas( 64 ) const std::array<int, 16> s_zeroTailMask =
+{
+	-1,-1,-1,-1,-1,-1,-1,-1,
+	0, 0, 0, 0, 0, 0, 0, 0,
+};
+
+namespace
+{
+	__forceinline float horizontalSum( __m256 vec )
+	{
+		__m128 v = _mm256_extractf128_ps( vec, 1 );
+		v = _mm_add_ps( v, _mm256_castps256_ps128( vec ) );
+		v = _mm_add_ps( v, _mm_movehl_ps( v, v ) );
+		v = _mm_add_ss( v, _mm_movehdup_ps( v ) );
+		return _mm_cvtss_f32( v );
+	}
+}
+
+void norm( float* rdi, float* temp, const float* rsi, size_t length )
+{
+	assert( (size_t)temp % 32 == 0 );
+	const float* rsiEndAligned = rsi + ( length & maskAlign8 );
+	const size_t rem = length % 8;
+
+	// First pass: copy to temp buffer, and compute the sum; computeVectorSum() in HLSL
+	__m256 sum = _mm256_setzero_ps();
+	float* t;
+	for( t = temp; rsi < rsiEndAligned; rsi += 8, t += 8 )
+	{
+		__m256 v = _mm256_loadu_ps( rsi );
+		sum = _mm256_add_ps( sum, v );
+		_mm256_store_ps( t, v );
+	}
+	float* const tEndAligned = t;
+	if( 0 != rem )
+	{
+		__m256 v = loadPartial( rsi, rem );
+		sum = _mm256_add_ps( sum, v );
+		_mm256_store_ps( t, v );
+		t += 8;
+	}
+
+	const float lengthFloat = (float)(int)length;
+	const float meanScalar = horizontalSum( sum ) / lengthFloat;
+	const __m256 mean = _mm256_set1_ps( meanScalar );
+
+	// Second pass, offsetAndComputeSumSquares() in HLSL
+	sum = _mm256_setzero_ps();
+	for( t = temp; t < tEndAligned; t += 8 )
+	{
+		__m256 v = _mm256_load_ps( t );
+		v = _mm256_sub_ps( v, mean );
+		_mm256_store_ps( t, v );
+		sum = _mm256_fmadd_ps( v, v, sum );
+	}
+	if( 0 != rem )
+	{
+		__m256 v = _mm256_load_ps( t );
+		v = _mm256_sub_ps( v, mean );
+		v = _mm256_and_ps( v, loadTailMaskFloats( rem ) );
+		_mm256_store_ps( t, v );
+		sum = _mm256_fmadd_ps( v, v, sum );
+	}
+
+	// Final pass: scale, and copy from temporary buffer into the destination row
+
+	constexpr float eps = 1e-5f; // TODO: make this a parameter
+	const float scaleScalar = 1.0f / std::sqrtf( horizontalSum( sum ) / lengthFloat + eps );
+	const __m256 scale = _mm256_set1_ps( scaleScalar );
+
+	for( t = temp; t < tEndAligned; t += 8, rdi += 8 )
+	{
+		__m256 v = _mm256_load_ps( t );
+		v = _mm256_mul_ps( v, scale );
+		_mm256_storeu_ps( rdi, v );
+	}
+	if( 0 != rem )
+	{
+		__m256 v = _mm256_load_ps( t );
+		v = _mm256_mul_ps( v, scale );
+		storePartial( rdi, v, rem );
+	}
+}
+
+void fmaRepeatRow( float* rdi, size_t len, const float* w, const float* b, size_t lenPattern )
+{
+	float* rdiEndAligned = rdi + ( len & maskAlign8 );
+	const size_t rem = len % 8;
+
+	if( 1 == lenPattern )
+	{
+		const __m256 v1 = _mm256_broadcast_ss( w );
+		const __m256 v2 = _mm256_broadcast_ss( b );
+		for( ; rdi < rdiEndAligned; rdi += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			v = _mm256_fmadd_ps( v, v1, v2 );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			v = _mm256_fmadd_ps( v, v1, v2 );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+	}
+	else if( len == lenPattern )
+	{
+		for( ; rdi < rdiEndAligned; rdi += 8, w += 8, b += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			__m256 v1 = _mm256_loadu_ps( w );
+			__m256 v2 = _mm256_loadu_ps( b );
+			v = _mm256_fmadd_ps( v, v1, v2 );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			__m256 v1 = _mm256_maskload_ps( w, mask );
+			__m256 v2 = _mm256_maskload_ps( b, mask );
+			v = _mm256_fmadd_ps( v, v1, v2 );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+	}
+	else
+	{
+		// TODO: implement if this actually happens
+		throw E_NOTIMPL;
+	}
+}
+
+void __vectorcall addRepeatScaleRow( float* rdi, size_t len, const float* b, size_t lenPattern, const __m256 scale )
+{
+	float* rdiEndAligned = rdi + ( len & maskAlign8 );
+	const size_t rem = len % 8;
+
+	if( 1 == lenPattern )
+	{
+		const __m256 v2 = _mm256_broadcast_ss( b );
+		for( ; rdi < rdiEndAligned; rdi += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			v = _mm256_add_ps( v, v2 );
+			v = _mm256_mul_ps( v, scale );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			v = _mm256_add_ps( v, v2 );
+			v = _mm256_mul_ps( v, scale );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else if( len == lenPattern )
+	{
+		for( ; rdi < rdiEndAligned; rdi += 8, b += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			__m256 v2 = _mm256_loadu_ps( b );
+			v = _mm256_add_ps( v, v2 );
+			v = _mm256_mul_ps( v, scale );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			__m256 v2 = _mm256_maskload_ps( b, mask );
+			v = _mm256_add_ps( v, v2 );
+			v = _mm256_mul_ps( v, scale );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else
+	{
+		// TODO: implement if this actually happens
+		throw E_NOTIMPL;
+	}
+}
+
+void addRepeatRow( float* rdi, size_t len, const float* b, size_t lenPattern )
+{
+	float* rdiEndAligned = rdi + ( len & maskAlign8 );
+	const size_t rem = len % 8;
+
+	if( 1 == lenPattern )
+	{
+		const __m256 v2 = _mm256_broadcast_ss( b );
+		for( ; rdi < rdiEndAligned; rdi += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			v = _mm256_add_ps( v, v2 );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			v = _mm256_add_ps( v, v2 );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else if( len == lenPattern )
+	{
+		for( ; rdi < rdiEndAligned; rdi += 8, b += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			__m256 v2 = _mm256_loadu_ps( b );
+			v = _mm256_add_ps( v, v2 );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			__m256 v2 = _mm256_maskload_ps( b, mask );
+			v = _mm256_add_ps( v, v2 );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else
+	{
+		// TODO: implement if this actually happens
+		throw E_NOTIMPL;
+	}
+}
+
+namespace
+{
+	__forceinline __m256 gelu( __m256 x, const DirectCompute::LookupTablesData& lookup )
+	{
+		__m128i iv = _mm256_cvtps_ph( x, 0 );
+		alignas( 16 ) std::array<uint16_t, 8> arr;
+		_mm_store_si128( ( __m128i* )arr.data(), iv );
+		for( uint16_t& a : arr )
+			a = lookup.gelu[ a ];
+		iv = _mm_load_si128( ( __m128i* )arr.data() );
+		return _mm256_cvtph_ps( iv );
+	}
+}
+
+void addRepeatGeluRow( float* rdi, size_t len, const float* b, size_t lenPattern, const DirectCompute::LookupTablesData& lookup )
+{
+	float* rdiEndAligned = rdi + ( len & maskAlign8 );
+	const size_t rem = len % 8;
+
+	if( 1 == lenPattern )
+	{
+		const __m256 v2 = _mm256_broadcast_ss( b );
+		for( ; rdi < rdiEndAligned; rdi += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			v = _mm256_add_ps( v, v2 );
+			v = gelu( v, lookup );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			v = _mm256_add_ps( v, v2 );
+			v = gelu( v, lookup );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else if( len == lenPattern )
+	{
+		for( ; rdi < rdiEndAligned; rdi += 8, b += 8 )
+		{
+			__m256 v = _mm256_loadu_ps( rdi );
+			__m256 v2 = _mm256_loadu_ps( b );
+			v = _mm256_add_ps( v, v2 );
+			v = gelu( v, lookup );
+			_mm256_storeu_ps( rdi, v );
+		}
+		if( 0 != rem )
+		{
+			const __m256i mask = loadTailMaskInt( rem );
+			__m256 v = _mm256_maskload_ps( rdi, mask );
+			__m256 v2 = _mm256_maskload_ps( b, mask );
+			v = _mm256_add_ps( v, v2 );
+			v = gelu( v, lookup );
+			_mm256_maskstore_ps( rdi, mask, v );
+		}
+		return;
+	}
+	else
+	{
+		// TODO: implement if this actually happens
+		throw E_NOTIMPL;
+	}
+}
+
+void __vectorcall scaleRow( float* rdi, size_t len, const __m256 scale )
+{
+	float* rdiEndAligned = rdi + ( len & maskAlign8 );
+	const size_t rem = len % 8;
+	for( ; rdi < rdiEndAligned; rdi += 8 )
+	{
+		__m256 v = _mm256_loadu_ps( rdi );
+		v = _mm256_mul_ps( v, scale );
+		_mm256_storeu_ps( rdi, v );
+	}
+	if( 0 != rem )
+	{
+		const __m256i mask = loadTailMaskInt( rem );
+		__m256 v = _mm256_maskload_ps( rdi, mask );
+		v = _mm256_mul_ps( v, scale );
+		_mm256_maskstore_ps( rdi, mask, v );
+	}
+}
+
+namespace
+{
+	using DirectCompute::LookupTablesData;
+
+	__forceinline float horizontalMax( __m256 vec )
+	{
+		__m128 v = _mm256_extractf128_ps( vec, 1 );
+		v = _mm_max_ps( v, _mm256_castps256_ps128( vec ) );
+		v = _mm_max_ps( v, _mm_movehl_ps( v, v ) );
+		v = _mm_max_ss( v, _mm_movehdup_ps( v ) );
+		return _mm_cvtss_f32( v );
+	}
+
+	__forceinline float _cvtsh_ss( uint16_t f16 )
+	{
+		__m128i i = _mm_cvtsi32_si128( f16 );
+		__m128 f = _mm_cvtph_ps( i );
+		return _mm_cvtss_f32( f );
+	}
+
+	__forceinline uint16_t _cvtss_sh( float f, int rounding )
+	{
+		assert( 0 == rounding );
+		__m128 v = _mm_set_ss( f );
+		__m128i i = _mm_cvtps_ph( v, 0 );
+		return (uint16_t)(uint32_t)_mm_cvtsi128_si32( i );
+	}
+}
+
+const LookupTablesData& getLookupTables()
+{
+	static const std::unique_ptr<LookupTablesData> res = std::make_unique<LookupTablesData>();
+	return *res;
+}
+
+void softMax( float* rdi, size_t length, const float inputScale )
+{
+	float* const rdiBegin = rdi;
+	float* const rdiEndAligned = rdi + ( length & maskAlign8 );
+	const size_t remainder = length % 8;
+	// First pass, compute maximum
+	__m256 max = _mm256_set1_ps( -INFINITY );
+	for( rdi = rdiBegin; rdi < rdiEndAligned; rdi += 8 )
+	{
+		__m256 v = _mm256_loadu_ps( rdi );
+		max = _mm256_max_ps( max, v );
+	}
+	__m256i tailMask;
+	if( 0 != remainder )
+	{
+		tailMask = loadTailMaskInt( remainder );
+		__m256 v = _mm256_maskload_ps( rdi, tailMask );
+		v = _mm256_max_ps( max, v );
+		max = _mm256_blendv_ps( max, v, _mm256_castsi256_ps( tailMask ) );
+	}
+
+	// Second pass: apply initial scale, compute the exponent, and compute total sum over the row
+	const LookupTablesData& lookup = getLookupTables();
+	const float maxScalar = horizontalMax( max );
+
+	float* const rdiEnd = rdiBegin + length;
+	double sum = 0;
+	for( rdi = rdiBegin; rdi < rdiEnd; rdi++ )
+	{
+		// Possible to vectorize, but relatively hard
+		// An easy way is upcast the complete lookup table to FP32 and then use two _mm256_i32gather_ps instructions per iteration
+		// However, that instruction is from AVX2 set. Let's hope this loop won't be a bottleneck.
+		float f = *rdi;
+		if( f != -INFINITY )
+		{
+			f = ( f - maxScalar ) * inputScale;
+			uint16_t f16 = _cvtss_sh( f, 0 );
+			f16 = lookup.exponent[ f16 ];
+			f = _cvtsh_ss( f16 );
+			sum += f;
+		}
+		else
+			f = 0;
+
+		*rdi = f;
+	}
+
+	// Final pass: apply the final scale
+	const __m256 finalScale = _mm256_set1_ps( (float)( 1.0 / sum ) );
+	for( rdi = rdiBegin; rdi < rdiEndAligned; rdi += 8 )
+	{
+		__m256 v = _mm256_loadu_ps( rdi );
+		v = _mm256_mul_ps( v, finalScale );
+		_mm256_storeu_ps( rdi, v );
+	}
+	if( 0 != remainder )
+	{
+		__m256 v = _mm256_maskload_ps( rdi, tailMask );
+		v = _mm256_mul_ps( v, finalScale );
+		_mm256_maskstore_ps( rdi, tailMask, v );
+	}
+}
+
+void floatsUpcast( float* rdi, const uint16_t* rsi, size_t length )
+{
+	const uint16_t* rsiEndAligned = rsi + ( length & maskAlign8 );
+	const size_t rem = length % 8;
+
+	for( ; rsi < rsiEndAligned; rsi += 8, rdi += 8 )
+		_mm256_storeu_ps( rdi, load8( rsi ) );
+
+	if( 0 != rem )
+	{
+		__m256 v = loadPartial( rsi, rem );
+		_mm256_maskstore_ps( rdi, loadTailMaskInt( rem ), v );
+	}
+}
+
+void floatsDowncast( uint16_t* rdi, const float* rsi, size_t length )
+{
+	const float* rsiEndAligned = rsi + ( length & maskAlign8 );
+	size_t rem = length % 8;
+
+	for( ; rsi < rsiEndAligned; rsi += 8, rdi += 8 )
+	{
+		__m256 vf = _mm256_loadu_ps( rsi );
+		__m128i vi = _mm256_cvtps_ph( vf, 0 );
+		store16( rdi, vi );
+	}
+
+	if( 0 != rem )
+	{
+		__m256 vf = _mm256_maskload_ps( rsi, loadTailMaskInt( rem ) );
+		__m128i vi = _mm256_cvtps_ph( vf, 0 );
+		for( size_t i = 0; i < rem; i++, rdi++ )
+		{
+			*rdi = (uint16_t)(uint32_t)_mm_cvtsi128_si32( vi );
+			vi = _mm_srli_si128( vi, 2 );
+		}
+	}
+}
+
+void addRowInPlace( float* rdi, const float* rsi, size_t length )
+{
+	const float* rdiEndAligned = rdi + ( length & maskAlign8 );
+	size_t rem = length % 8;
+
+	for( ; rdi < rdiEndAligned; rdi += 8, rsi += 8 )
+	{
+		__m256 a = _mm256_loadu_ps( rdi );
+		__m256 b = _mm256_loadu_ps( rsi );
+		a = _mm256_add_ps( a, b );
+		_mm256_storeu_ps( rdi, a );
+	}
+
+	if( 0 != rem )
+	{
+		const __m256i mask = loadTailMaskInt( rem );
+		__m256 a = _mm256_maskload_ps( rdi, mask );
+		__m256 b = _mm256_maskload_ps( rsi, mask );
+		a = _mm256_add_ps( a, b );
+		_mm256_maskstore_ps( rdi, mask, a );
+	}
+}
+
+void addRow( float* rdi, const float* a, const float* b, size_t length )
+{
+	const float* aEndAligned = a + ( length & maskAlign8 );
+	size_t rem = length % 8;
+
+	for( ; a < aEndAligned; a += 8, b += 8, rdi += 8 )
+	{
+		__m256 x = _mm256_loadu_ps( a );
+		__m256 y = _mm256_loadu_ps( b );
+		x = _mm256_add_ps( x, y );
+		_mm256_storeu_ps( rdi, x );
+	}
+
+	if( 0 != rem )
+	{
+		const __m256i mask = loadTailMaskInt( rem );
+		__m256 x = _mm256_maskload_ps( a, mask );
+		__m256 y = _mm256_maskload_ps( b, mask );
+		x = _mm256_add_ps( x, y );
+		_mm256_maskstore_ps( rdi, mask, x );
+	}
+}
\ No newline at end of file
diff --git a/Whisper/CPU/simdUtils.h b/Whisper/CPU/simdUtils.h
new file mode 100644
index 0000000..a7a4bac
--- /dev/null
+++ b/Whisper/CPU/simdUtils.h
@@ -0,0 +1,82 @@
+#pragma once
+#include <immintrin.h>
+
+void addF16to32( float* rdi, const uint16_t* a, const uint16_t* b, size_t length );
+void addF16to32( float* rdi, const uint16_t* a, const float* b, size_t length );
+
+class AlignedSpan
+{
+	float* pointer;
+
+public:
+	AlignedSpan( void* data )
+	{
+		size_t i = (size_t)data;
+		constexpr size_t mask32 = ~(size_t)31;
+		i = ( i + 31 ) & mask32;
+		pointer = (float*)i;
+	}
+
+	operator float* ( ) { return pointer; }
+};
+
+inline size_t tempBufferForFloats( size_t count )
+{
+	// Round up by 8 to be able to use full-vector loads and stores
+	constexpr size_t mask8 = ~(size_t)7;
+	count = ( count + 7 ) & mask8;
+
+	// Add 32 more bytes to align the temporary buffer
+	return ( count * 4 ) + 32;
+}
+
+#define ALIGNED_SPAN( name, countFloats ) AlignedSpan name{ _alloca( tempBufferForFloats( countFloats ) ) }
+
+void norm( float* rdi, float* temp, const float* rsi, size_t length );
+
+void fmaRepeatRow( float* rdi, size_t len, const float* w, const float* b, size_t lenPattern );
+void __vectorcall addRepeatScaleRow( float* rdi, size_t len, const float* b, size_t lenPattern, const __m256 scale );
+void addRepeatRow( float* rdi, size_t len, const float* b, size_t lenPattern );
+void __vectorcall scaleRow( float* rdi, size_t len, const __m256 scale );
+
+namespace DirectCompute
+{
+	struct LookupTablesData;
+}
+const DirectCompute::LookupTablesData& getLookupTables();
+void addRepeatGeluRow( float* rdi, size_t len, const float* b, size_t lenPattern, const DirectCompute::LookupTablesData& lookup );
+
+void softMax( float* rdi, size_t length, const float inputScale );
+
+// A cache line-aligned array where first 8 elements have all bits set, last 8 elements are zeros
+extern const std::array<int, 16> s_zeroTailMask;
+
+// Load a tail mask as FP32 vector, for use with _mm256_and_ps or _mm256_blendv_ps instructions
+__forceinline __m256 loadTailMaskFloats( size_t remainder )
+{
+	assert( remainder > 0 && remainder < 8 );
+	const float* rsi = (const float*)&s_zeroTailMask;
+	rsi += 8;
+	return _mm256_loadu_ps( rsi - remainder );
+}
+
+// Load a tail mask as int32 vector, for use with _mm256_maskstore_ps instruction
+template<bool assertIncomplete = true>
+__forceinline __m256i loadTailMaskInt( size_t remainder )
+{
+	if constexpr( assertIncomplete )
+		assert( remainder > 0 && remainder < 8 );
+	else
+		assert( remainder >= 0 && remainder <= 8 );
+
+	const int* rsi = (const int*)&s_zeroTailMask;
+	rsi += 8;
+	return _mm256_loadu_si256( ( const __m256i* )( rsi - remainder ) );
+}
+
+void floatsUpcast( float* rdi, const uint16_t* rsi, size_t length );
+
+void floatsDowncast( uint16_t* rdi, const float* rsi, size_t length );
+
+void addRowInPlace( float* rdi, const float* rsi, size_t length );
+void addRow( float* rdi, const float* a, const float* b, size_t length );
\ No newline at end of file