Using signed integer loop counters is leading to some extra work outside the loop, but it doesn't seem to be hurting the inner-most loop in the compiler's asm output. Still, you might as well use size_t for memory sizes; that's what it's for and usually compiles efficiently. OTOH, limiting your args to signed makes it easier to express the condition that the loop should run zero times (because size - unroll width can be negative, instead of having to code it like asm, with an explicit compare and then a do{}while() loop.)
numColsPad should be called rowStride or stride. numColsPad is ambiguous whether it means number of columns including padding, or whether it's the width of the padding! Stride is a well-established term for the distance in memory between rows. stride - cols doesn't even have to be small: separating the storage layout from the width lets you pass a sub-image to a function without copying the data to another buffer.
numColsSse is an awkward variable name. SSE should be all-caps. Or better, just call it something else, like colVecs. I'm not a fan of your other variable names, or camelCase in general. HorizontalMaxSse is also a silly name. Max of what? Call it HorizontalMax_ps to indicate that it's packed-single-precision floats.
You should use const float* for the read-only input array. You could also return min by value and only use the output-pointer for max. It would be better C++ style to declare your variables as you use them, instead of C89-style at the top of the function.
The two biggest performance things I see are:
Handling the tail of every row with a 2nd loop over all the rows is horrible (for(jj = numColsSse + 1; ...). That data will be hot in L1D cache from HW prefetch while doing a row, so you should grab it then instead of waiting until later when it's probably evicted again. Striding over the rows to load a couple columns is a bad access pattern which will hurt a lot for tall and narrow image.
You also don't need to do it scalar. You can accumulate tailMin and tailMax vectors, and after the fact take only the elements that are actually part of the image. (e.g. loop over 1 to 4 elements of the vector, depending on numCols%4. Not 0-3, because that would mean we'd want to do it conditionally to avoid wasting work.)
If you need to avoid loading past the end of the last row, then peel the last row out of the main loop and do it specially.
It makes the cleanup trivial if we use an unaligned potentially-overlapping load that gets the last 4 elements of a row (then we don't have to ignore any elements of separate tail accumulators). Overlap is fine because x=min(x,y) is idempotent. Unlike a sum reduction, looking at the same element multiple times doesn't change the result. If we're clever, we can use this to simplify the unrolling.
OTOH, branching based on numCols & 7 would probably be better to avoid doing a fully-overlapped vector, since we bottleneck on min/maxps throughput, not on the front-end. Still, this is pretty good and makes the source simple. Optimize further however you like.
Unroll the inner vector loop with multiple accumulators. maxps has 3c or 4c latency, 1c or 0.5c throughput on recent Intel CPUs (http://agner.org/optimize/, and see also this SO Q&A about using multiple accumulators for FMA.)
On Skylake, you need 8 accumulators to keep 8 min/maxps instructions in flight at once to max out execution-unit throughput (it runs FP add/cmp/min/max on the FMA units on p0/p1, same as mul and fma). You only have 2, one each for min and max. On Broadwell (3c latency, 1c throughput, runs on port1 only) this will in theory get you 2/3rds of theoretical max throughput.
This is what MSVC CL19 does for the inner loop:
$LL7@ExtractMin:
movups xmm0, XMMWORD PTR [rdx]
lea rdx, QWORD PTR [rdx+16] ; ADD would be better, but this probably doesn't create resource conflicts
minps xmm1, xmm0
maxps xmm2, xmm0
sub r8, 1 ; will macro-fuse with JNE on Intel SnB-family
jne SHORT $LL7@ExtractMin
That's 5 fused-domain uops, so it can issue at about 1 iteration per 1.25 clocks. So actually even on Haswell/Broadwell, min/maxps latency is a bigger bottleneck than CPU front-end throughput (4 uops per clock max, or slightly less in loops that aren't a multiple of 4 uops in some cases.) Still, unrolling will let out-of-order execution see the loads earlier.
OTOH, unrolling introduces more code to get the iteration count right (especially if you need your code to work for very narrow images where you can't assume the width is at least one full vector).
This can be outside the outer loop if the whole nested loop is duplicated. (Either in the source or in the compiler output if it decides to factor out a check). Or we can do something simple and branchless like always doing a first vector and then starting j from 0 or 4 depending on whether numCols is a multiple of SSE_STRIDE*2 or not (ignoring the low 2 bits). Overlap is fine because min and max are idempotent. (I haven't thought through the possibilities with NaN; _mm_min_ps() and max are not commutative without -ffast-math or equivalent for other compilers. You are using them with the operand order that allows the compiler to avoid extra movaps instructions.)
Here's an example of most of those ideas. I'm only unrolling by 2, but you need to unroll by at least 4 to max out performance on Skylake (for cases where you don't bottleneck on memory).
source + asm on Godbolt. It looks good for MSVC and gcc and clang. Clang does further unrolling of the inner loop, but without using more accumulators so it's useless for speeding up Skylake.
Returning a 2-member struct by value ends up taking slightly more code than just storing twice, so it's probably not a win. I'd recommend picking one method (or just returning min by value as a float and max by pointer), instead of storing both and packing both into a struct. The caller will have to unpack the struct into separate registers in most cases, so it's not that great.
#include <immintrin.h>
struct minmax { float min, max; };
// pre-condition: numCols >= 8
// numRows >= 1
minmax ExtractMinMax(float* vMinMax, const float* mI, int numRows, int numCols, int rowStride)
{
static const int VEC_WIDTH = 4;
//if (numCols < 2*VEC_WIDTH) {
// return alternate_version(vMinMax, mI, numRows, numCols, rowStride);
//}
// if you want to support numCols < 8, a simple way is _mm_load1_ps(mI)
// but you'd also have to change the non-multiple-of-8 handling in the loop
__m128 min0 = _mm_loadu_ps(mI);
__m128 min1 = _mm_loadu_ps(mI + VEC_WIDTH);
__m128 max0 = min0, max1 = min1;
// TODO: skip the first 2 vectors of the first row
int i=0;
do { // let the compile know this loop runs at least once
// manually hoisting the row pointer mostly helps human readability; compilers usually already do this
const float *row = &mI[i * rowStride];
for (int j=0 ; j < numCols - VEC_WIDTH*2; j += VEC_WIDTH*2) {
// for large unroll counts, a function or macro is useful. (not recommended, don't go nuts)
// for small unrolls, manual repetition is usually the most readable.
__m128 pxI0 = _mm_loadu_ps(row+j);
min0 = _mm_min_ps(min0, pxI0);
max0 = _mm_max_ps(max0, pxI0);
__m128 pxI1 = _mm_loadu_ps(row+j + VEC_WIDTH);
min1 = _mm_min_ps(min1, pxI1);
max1 = _mm_max_ps(max1, pxI1);
}
// loop stops with 1 to 2*VEC_WIDTH elements left to do, depending on numCols & (2*VEC_WIDTH - 1)
// last 2 vectors ending at the end of the row, could overlap with above
// but won't for multiple-of-8 numCols
__m128 tail = _mm_loadu_ps(row + numCols - VEC_WIDTH*2);
min0 = _mm_min_ps(min0, tail);
max0 = _mm_max_ps(max0, tail);
tail = _mm_loadu_ps(row + numCols - VEC_WIDTH);
min1 = _mm_min_ps(min1, tail);
max1 = _mm_max_ps(max1, tail);
// This isn't really optimal; 2 unaligned loads is more chances to CL-split or 4k-split
// assuming the start of each row is aligned
// but it does avoid any branching on the low bits of numCols.
} while(++i < numRows);
// reduce the multiple accumulators to one
min0 = _mm_min_ps(min0, min1);
max0 = _mm_max_ps(max0, max1);
// horizontal min/max.
float min = HorizontalMin_ps(min0);
float max = HorizontalMax_ps(max0);
vMinMax[0] = min;
vMinMax[1] = max;
return { min, max };
}
Using AVX if available could help, too, but probably not if you're totally memory bottlenecked from looping over giant images doing not much work for the amount of memory traffic. Still, AVX can help saturate memory bandwidth better because the same number of in-flight loads covers twice as much data. (But requests from L1 to L2, and L2 to L3, are tracked in whole cache lines.)
Try to combine the min/max calculation into whatever wrote the image in the first place, to increase computational intensity. Doing a separate pass for min/max isn't great. (Although it's twice as good as doing separate passes over your data for min and for max.)
Using multiple threads can help even if you're memory-bound, especially on Intel CPUs with more than quad cores. IDK how bad single-threaded memory performance is on your 6-core Broadwell i7-6800k, but on many-core Xeons the extra latency to L3 and DRAM hurts single-threaded bandwidth a lot. (bandwidth = max concurrency / latency). See the Latency-bound platforms section of this answer for more details.
My quad-core Skylake (i7-6700k) can use most of the total memory bandwidth with only one thread, but that's with a quite high core / uncore clock speed.
Also, with enough cores for a medium-size image, it might start fitting in L2 cache (since each core only looks at part of it). But 2k*2k * 4B is 15.25MiB, which is larger than the 15MB L3 cache on your CPU.
If you're going to use this with unaligned input pointers very often (e.g. a sub-image), it might make sense to do the unaligned head of each row separately, and a potentially-overlapping aligned loop over the main part of each row. cache-line splits won't be a problem on recent Intel (you're only going to do 1 load per clock max), but 4k-splits will still cause bubbles on pre-Skylake. Like 34 cycle load-use latency on Haswell for a 4k-split instead of only 11 for a cache-line split.
NaN behaviour: this will ignore NaN in the input (unless it's in the initial values that we start with.
max and _mm_min_ps() are not commutative without -ffast-math (or equivalent for other compilers). You are using them with the operand order that allows the compiler to avoid extra movaps instructions, so that's a good thing.
HorizontalMinSse(and Max), andSSE_STRIDE? (The horizontal code probably won't matter on any normal image size, though, since it's only run once per image) \$\endgroup\$numCols % SSE_STRIDEis non-zero, and those will all miss in cache because the vector loop over the main part of each row touched so much memory. Anyway, you should look at your CPU load to figure out if MATLAB is multi-threading. \$\endgroup\$