Open jungjanos opened 4 years ago
I´ll have a look later, thanks. FYI, my code runs without errors here (both scalar and vectorized).
Hey,
I have also implemented an AVX version: A accumulate the iteracion counts for the currently processed eight pixels in a Vector256
public static unsafe int[] VectorMandel_Avx_256(int resolutionY, int maxIterations)
{
const int VECTOR_SIZE = 8;
// ensure that resolutionX is divisible by SIMD-width
// ensure that resolutionY / 2.0 * 3.5 is integer
// this simplifies calculations
if (resolutionY % (4 * VECTOR_SIZE) != 0)
throw new ArgumentException($"resolutionY must be divisible by {VECTOR_SIZE * 4}");
int resolutionX = resolutionY / 4 * 7;
// distance between neighbouring points on the X-Y plane
float step = 2f / resolutionY;
// stepper vectors to easily step in Y and X directions
var yStepVec = Vector256.Create(-step);
var xStepVec = Vector256.Create(VECTOR_SIZE * step);
// always contains the current Y-axis position. We begin at the top most horizontal line
// will be updated inside loop. All SIMD-components of currentYvec are always the same
var currentYvec = Vector256.Create(1f);
var result = new int[resolutionX * resolutionY];
fixed (int* resultFixedPtr = result)
{
int* resultPtr = resultFixedPtr;
// outer loop: iterate along Y axis from top to bottom, one step at a time
for (int i = 0; i < resolutionY; i++)
{
// initialize the starting X vector to the leftmost X positions
var currentXvec = Vector256.Create(
-2.5f, -2.5f + step, -2.5f + 2 * step, -2.5f + 3 * step,
-2.5f + 4 * step, -2.5f + 5 * step, -2.5f + 6 * step, -2.5f + 7 * step
);
// iterate along the X axis with "VECTOR_SIZE" number of elements at a time
// the value of "currentXvec" is updated at the end of the loop to always represent the working elements
for (int j = 0; j <= resolutionX - VECTOR_SIZE; j += VECTOR_SIZE)
{
int iteration = 0;
var iterationCounter = Vector256.Create(0f);
Vector256<float> xVec, yVec;
var xVecSquare = Vector256.Create(0f);
var yVecSquare = Vector256.Create(0f);
var sumVecSquare = Vector256.Create(0f);
var fourVec = Vector256.Create(4f);
var oneVec = Vector256.Create(1f);
while (iteration < maxIterations)
{
var test = Avx.Compare(Avx.Add(xVecSquare, yVecSquare), fourVec,FloatComparisonMode.OrderedLessThanOrEqualNonSignaling);
// no components inside tolerance, no work left to do
if (Avx.MoveMask(test) == 0)
break;
// mask the incrementer vector to increment only the positions still inside Mandelbrot-tolerance
var selectiveAdd = Avx.And(oneVec, test);
iterationCounter = Avx.Add(iterationCounter, selectiveAdd);
iteration++;
// set up values for next iteration
xVec = Avx.Add(Avx.Subtract(xVecSquare, yVecSquare), currentXvec);
yVec = Avx.Add(Avx.Subtract(sumVecSquare, Avx.Add(xVecSquare, yVecSquare)), currentYvec);
xVecSquare = Avx.Multiply(xVec, xVec);
yVecSquare = Avx.Multiply(yVec, yVec);
sumVecSquare = Avx.Multiply(Avx.Add(xVec, yVec), Avx.Add(xVec, yVec));
}
// at this point the iterations counts for a whole vector has been calculated
// resulting values are stored
Avx.Store(resultPtr, Avx.ConvertToVector256Int32WithTruncation(iterationCounter));
resultPtr += VECTOR_SIZE;
// move along the X axis by VECTOR_SIZE step-s
currentXvec = Avx.Add(currentXvec, xStepVec);
}
// move along the Y axis by one step
currentYvec = Avx.Add(currentYvec, yStepVec);
}
}
return result;
}
I had a look ate the SSE code and run it: if I do a pixel-by-bixel comparison (compare the result arrays, not the bitmaps, too lazy to get the nuget package... ) I get differences, I suspect due to floating point deltas... not an issue with Mandelbrot but could be an issue for production code. That´s one of the reasons I pre-calculate the x and y values since I test the benchmarks by comparing the scalar and vector outputs. Did not look ate the Avx code yet. I´ll see if I can eliminate the differences and include your code and benchmarks if that´s OK with you. Cheers!
You are right, my scalar and my vector codes do not return the same results. (the SSE and AVX versions do get the same results though). I have checked now, the differences are relativelly few, but definitively are there.
I have debuged the code, I also have identified the reason: The reason is, that (in most cases) floating point numbers are only an approximation of the mathematical value (e.g. 1.1 decimal can not be represented with absolute accuracy in float) thus the results of floating point arithmetic are only an approximation of the real results.
The important thing is, that when you do floating point arithmetic e.g.: (a-b-c+d) than the execution order of the arithmetic matters!!!. a-b-c+d is not always equal to (a-(b+c)+d).
Exactly that is the difference between my scalar and vector codes: if you substitute the lines:
// yVec = Sse42.Add(Sse42.Subtract(sumVecSquare, Sse42.Add(xVecSquare, yVecSquare)), currentYvec); <---- order not conforming to scalar method
yVec = Sse42.Add(Sse42.Subtract(Sse42.Subtract(sumVecSquare, xVecSquare), yVecSquare), currentYvec); <---- order conforming to scalar method
than the results will be 100% identical to the scalar code results. (I have checked it)
Hi!
We have had contact on Reddit. I am pretty sure that your SIMD code doesnt work correctly and it also has design fault.
I have rewritten both the Scalar and the SIMD code, mainly to make it simpler. Take what you think useful.
I very much like your buffer size to x-y-dimensions calculations but I think it is an unnecessary overkill. So I just take the Y-resolution (number of vertical pixels) and calculate the X-resolution (horizontal pixels). Instead of resolution correction, I only accept y-resolution which results in whole-number X-resolution.
SIMD: Where I think your implementation is problematic are the following items:
You have in-memory arrays for x and y values which are read into Vectors. That is unnecessary and thus suboptimal as SIMD code can easily get memory-bound. You can calculate the running X and Y vectors purely from registers (setting initial value and increment with step size) without touching memory.
I also get at
testSpan[resVectorNumber] = Avx.Add(xSquVec, ySquVec);
an index out of range exceptionI have used only SSE4.2 as my machine is not AVX2 capable. If you add
SixLabors.ImageSharp;
from NuGet you can actually get an image of the calculated fractal!