On performance of matrix transposition using omatcopy()

May be somebody would like to improve the performance of omatcopy(), so here's some info:

Forget to tell about it for more than a year.... One may write a code to transpose a matrix in a "sequential reading" manner as well as in a "sequential writing" manner. I tried and timed both of these approaches and compared them to the OpenBLAS's omatcopy() function. Here's the typical results on my AMD PhenomII X6 1090T hardware (one may try it on their own hardware if compiles the code below or dare to build nntl's tests project and run it with --gtest_filter=*.mTranspose parameter):

******* checking mTranspose() variations over 100x100 matrix (10000 elements) **************
sread:     8.521 us
swrite:    9.119 us
OBLAS:     8.906 us
******* checking mTranspose() variations over 1000x10 matrix (10000 elements) **************
sread:     8.569 us
swrite:    9.401 us
OBLAS:     9.119 us
******* checking mTranspose() variations over 10x1000 matrix (10000 elements) **************
sread:     8.825 us
swrite:    9.228 us
OBLAS:     8.469 us
******* checking mTranspose() variations over 1000x100 matrix (100000 elements) **************
sread:   158.612 us
swrite:  166.583 us
OBLAS:   126.275 us
******* checking mTranspose() variations over 100x1000 matrix (100000 elements) **************
sread:   166.476 us
swrite:  135.253 us           <----------------
OBLAS:   160.407 us
******* checking mTranspose() variations over 10000x10 matrix (100000 elements) **************
sread:   250.123 us
swrite:   96.029 us           <----------------
OBLAS:   265.792 us
******* checking mTranspose() variations over 10x10000 matrix (100000 elements) **************
sread:   118.552 us
swrite:  201.491 us
OBLAS:    90.101 us
******* checking mTranspose() variations over 10000x100 matrix (1000000 elements) **************
sread:     7.519 ms
swrite:    2.308 ms           <----------------
OBLAS:     7.495 ms
******* checking mTranspose() variations over 100x10000 matrix (1000000 elements) **************
sread:     2.321 ms
swrite:    6.909 ms
OBLAS:     2.122 ms

Here's the code excerpts:

//corresponds to "sread" entry in timing table
void mTranspose_seq_read(const realmtx_t& src, realmtx_t& dest) noexcept {
	NNTL_ASSERT(src.rows() == dest.cols() && src.cols() == dest.rows());
	const ptrdiff_t sRows = src.rows(), sCols = src.cols();
	const auto dataCnt = src.numel();
	auto pSrc = src.data();
	const auto pSrcE = pSrc + dataCnt;
	auto pDest = dest.data();
	
	while (pSrc != pSrcE) {
		auto pD = pDest++;
		auto pS = pSrc;
		pSrc += sRows;
		const auto pSE = pSrc;
		while (pS != pSE) {
			*pD = *pS++;
			pD += sCols;
		}
	}
}

//corresponds to "swrite" entry in timing table
void mTranspose_seq_write(const realmtx_t& src, realmtx_t& dest) noexcept {
	NNTL_ASSERT(src.rows() == dest.cols() && src.cols() == dest.rows());
	const ptrdiff_t sRows = src.rows(), sCols = src.cols();
	const auto dataCnt = src.numel();
	auto pSrc = src.data();
	auto pDest = dest.data();
	const auto pDestE = pDest + dataCnt;

	while (pDest != pDestE) {
		auto pS = pSrc++;
		auto pD = pDest;		
		pDest += sCols;
		const auto pDE = pDest;
		while (pD != pDE) {
			*pD++ = *pS;
			pS += sRows;
		}
	}
}

//corresponds to OBLAS entry
void mTranspose_OpenBLAS(const realmtx_t& src, realmtx_t& dest) noexcept {
	const auto sRows = src.rows(), sCols = src.cols();
	math::b_OpenBLAS::omatcopy(true, sRows, sCols, real_t(1.0), src.data(), sRows, dest.data(), sCols);
}

math::b_OpenBLAS::omatcopy() function in that case is just a thunk to cblas_somatcopy()

template<typename sz_t, typename fl_t>
static typename std::enable_if_t< std::is_same< std::remove_pointer_t<fl_t>, float>::value >
omatcopy(const bool bTranspose, const sz_t rows, const sz_t cols, const fl_t alpha,
	const fl_t* pA, const sz_t lda, fl_t* pB, const sz_t ldb)
{
	cblas_somatcopy(CblasColMajor, bTranspose ? CblasTrans : CblasNoTrans,
		static_cast<blasint>(rows), static_cast<blasint>(cols), alpha,
		pA, static_cast<blasint>(lda), pB, static_cast<blasint>(ldb));
}

TL;DR: Looks like OpenBLAS uses only "sequential reading" method. It may be very beneficial to use the sequential writing approach on tall (rows cnt > cols cnt) matrices (CblasColMajor mode). Some processor's cache/memory line implementation details might also affect.

Start by reproducing the reported timings with nntl's tests and the --gtest_filter=*.mTranspose parameter, then inspect the cblas_somatcopy() entry point and its math::b_OpenBLAS::omatcopy() thunk. Compare sequential-reading and sequential-writing behavior for the matrix shapes shown; done means a measured performance improvement without changing transpose results.