Let's review this code purely from a performance angle, without a focus on style or anything else (in addition to optimization suggestions, Peter already mentioned several things in areas other than performance).
First, you can play with all the algorithm discussed here in this github repo. I compiled it on Linux but it should approximately work on Windows ...
The biggest performance problem (for small matrices), and a common mistake, is doing this:
_mm256_extract_epi32(vec_multi, 0) + _mm256_extract_epi32(vec_multi, 1) +_mm256_extract_epi32(vec_multi, 2) +_mm256_extract_epi32(vec_multi, 3) +_mm256_extract_epi32(vec_multi, 4) +_mm256_extract_epi32(vec_multi, 5) +_mm256_extract_epi32(vec_multi, 6) +...
One problem that immediately jumps out at me is the loop-carried dependency through addps, which has a latency of either 3 or 4 (depending on the processor) while there are not nearly enough instructions there to fill all that time, so it's lost throughput. The typical solution is unrolling and using multiple accumulators. There's too much stuff there for me ...
Taking Benefits of The Out-of-Order Execution Engine
You can also read about The Out-of-Order Execution Engine
in the "Intel® 64 and IA-32 Architectures Optimization Reference Manual"
http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-optimization-manual.pdf section the 2.1.2, and take benefits of it.
For example, ...
As suggested in the comments, you should profile the code to see where the slowdown is.
Algebra Is Your Friend
One small suggestion I can offer is some algebraic improvements.
A linear interpolation between a and b can be formulated like this:
result = a * interp + b * (1 - interp);
That's 2 adds and 2 multiplies. You can use some algebra ...
Without addressing performance concerns, some trivial observations:
#include <omp.h> is unnecessary. (You use OpenMP, but don't call any OpenMP functions.)
The return type of main() should be int, not void.
The code also compiles with clang (LLVM), if you omit the -masm=intel option.
zeromat() could simply be memset(C, 0, n * sizeof(double)).
The code seems to get more and more questionable as we read downward. Starting at the bottom:
inline auto generate_px8n_palette()
for (auto i = 0; i < 256; ++i)
Writing auto instead of explicitly std:...
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 ...
I think you're underestimating the importance of the code outside the inner loop, and the mis-predicts of exiting the loop with different numbers of counts. With that iteration-count profile, it's very unlikely that regular throughput bottlenecks like dep chains or ALU uops are the major factor in your total run-time.
You're almost certainly getting a lot ...
I put your code up on the Godbolt Compiler Explorer to check out the asm from the intrinsics version. You're right that even with gcc 5.3 or clang 3.8, there are spills / reloads in the inner loop. So you may actually get a speedup from hand-written asm here, if those store-forwarding round trips aren't hidden by out-of-order execution.
The problem with ...
For the time output in main():
Use std::clock_t as opposed to clock_t in C++.
You could have another variable, endTime, for the ending time.
The computed time could also have its own variable, elapsedTime, to be printed.
Only the clock time needs to be cast to a double, not the macro as well. You should also use static_cast<>() as this is C++.
array equ rdi ; pointer to array
size equ rsi ; array size (count of elements)
value equ xmm0 ; value to process with
I don't particularly care for these. It seems like you're trying to make the code look more like a high-level language, but it seems to me that it ends up ...
forget asm (for now at least) and stick with SSE intrinsics - you can concentrate on optimisation and let the compiler worry about the implementation details like register allocation, instruction scheduling and loop unrolling
use a decent compiler - e.g. Intel ICC typically generates better code than gcc in most cases, and Visual Studio ...
why not load 4 packed __m128 and then store 32 Pixels at once
I think that's 16 pixels, but it's a good plan. The pack instructions were used inefficiently (in the linked question that was not really an issue) and that would be improved by processing more elements simultaneously. For example, something like this:
#define SSE_STRIDE 16
save the registers we intend to alter, failure to do so causes problems when gcc -O3 is used
It should be possible to avoid the pushes and pops. If you are changing the values of the constraints, you cannot have them as just "inputs" (which is where this code currently has them). Quoting the docs:
Do not modify the contents of input-only operands (...
Gaussian blurs can also be split into a horizontal and vertical pass without loss of precision. It's actually one of the characteristics which makes them so popular. You may find that using several box-blur passes in each direction is actually slower than a single gaussian-blur pass in each direction (depending on how well you're able to optimise a gaussian ...
When I was composing this question, I kept thinking about what questions people were going to ask me: Did you think about A? Did you try B? Why didn't you go with the obvious thing and do C?
Normally, this is a good thing, since I can often answer my own questions without ever having to hit the "Send" button.
In this case, there were 2 questions buzzing ...
Something to look into: data access patterns. You go across all image lines, processing the first few pixels (were the kernel straddles the boundary), then again across all image lines, processing the bulk of the pixels, and then a third time across all image lines, processing the last few pixels. If you combine the outer loops into a single loop, walking ...
Get Rid of Unused Arguments
There are arguments passed to this function that are not used. It appears that the foreground array mF doesn't get read or modified, so remove it from the list of arguments as it's just confusing.
The names of your function's arguments are extremely cryptic. mO is the output buffer? If so, name it outputBuffer or ...
Small differences due to precision are expected and can usually be ignored.
12 shuffles like that are a bit much, though not necessarily avoidable, depending on whether there is AVX support. With AVX, it is better to literally broadcast from memory, rather than emulate broadcasting with a load and shuffles. Even though this means there will be more loads, ...
Thanks for all comments and tips! Especially @harold This is my final version with the main goal to eliminate cache misses by never loading columns into a vector.
It removes the costly set and load function by working with pointer loading and storing instead of element. By working row by row the cachemisses are few and memory blocks are used more fully.
There is no need to take the absolute value. If X is real, then X^2 is always positive.
Consider using hypot(x,y) rather than sqrt(x*x+y*y). It is simpler and more accurate.
double distanceSqrt(int x1, int x2, int y1, int y2)
double dx = x2 - x1;
double dy = y2 - y1;
return sqrt(dx*dx + dy*dy);
Your SSE code looks like a reasonably accurate translation of your original C, but that C code doesn't immediately jump out as the best way to do things. In particular, I don't see where you gain anything by taking the absolute value before you do your multiplication. Given that your're squaring the value immediately afterwards anyway, it would appear ...
One improvement would be to restrict the number of (random) memory accesses. First notice is that for each sourceY[fixed] there's also a memory read from sourceY[fixed+1];
I believe these should be at least combined to single 64-bit memory accesses.
A more crucial improvement would be at higher level: is the interpolation really random, or could the ...