Array Acces Speed

[SilverShaded]

Hello, i have a program which is heavy on array access, so i've been tyring to optimise the speed. This is important for the overall program which is highly iterative with lots of matrix manipulation. I was hoping to get any criticism of the code below or any further suggestions to speed it up more.

Current run times (in release mode) are respectively, in ms, the relative time differences are relatively constant, at the moment option 4 is looking like the best possible, which represent a significant increase in speed over option 2 which is the current method in the progam. Before i go changing a whole bunch of arrays are there any obvious improvements possible or better ways to do it or any problems with the code below that may give odd results?

elapsedMs1 128
elapsedMs2 29
elapsedMs3 228
elapsedMs4 7

 [TestMethod]
        [MethodImpl(MethodImplOptions.NoInlining)]
        public unsafe void TestDoubleDoubleArrayWithPts()
        {

            double[,] Array = new double[10000,1000];
            var watch = Stopwatch.StartNew();
            int UpperBound0 = Array.GetUpperBound(0);
            int UpperBound1 = Array.GetUpperBound(1);

            watch.Restart();
            for (int i = 0; i < Array.GetUpperBound(0); i++)
                for (int y = 0; y < Array.GetUpperBound(1); y++)
                    Array[i,y] = 10D * 10D;

            var elapsedMs1 = watch.ElapsedMilliseconds;

            watch.Restart();
            for (int i = 0; i < UpperBound0; i++)
                for (int y = 0; y < UpperBound1; y++)
                    Array[i, y] = 10D * 10D;

            var elapsedMs2 = watch.ElapsedMilliseconds;

            watch.Restart();
            fixed (double* ptr = Array)
            {
                for (int i = 0; i < Array.GetUpperBound(0); i++)
                    for (int y = 0; y < Array.GetUpperBound(1); y++)
                        *(ptr + i * Array.GetUpperBound(1) + y) = 10D * 10D;    
            }

            var elapsedMs3 = watch.ElapsedMilliseconds;

            watch.Restart();
            fixed (double* ptr = Array)
            {
                for (int i = 0; i < UpperBound0; i++)
                    for (int y = 0; y < UpperBound1; y++)
                        *(ptr + i*UpperBound1 + y) = 10D * 10D;      
            }

            var elapsedMs4 = watch.ElapsedMilliseconds;

            System.Windows.Forms.MessageBox.Show(elapsedMs1.ToString() + " " + elapsedMs2.ToString() + " " + elapsedMs3.ToString() + " " + elapsedMs4.ToString());
       
        }

 

[Skydiver]

I'm glad that is working out well for you.

Since BenchmarkDotNet is my shiny new toy, here's the results I got running with array size of 50x50 and 100x100:

BenchmarkDotNet=v0.12.1, OS=Windows 10.0.19041.746 (2004/?/20H1)
AMD Ryzen 7 PRO 2700U w/ Radeon Vega Mobile Gfx, 1 CPU, 8 logical and 4 physical cores
.NET Core SDK=5.0.102
  [Host]     : .NET Core 5.0.2 (CoreCLR 5.0.220.61120, CoreFX 5.0.220.61120), X64 RyuJIT
  DefaultJob : .NET Core 5.0.2 (CoreCLR 5.0.220.61120, CoreFX 5.0.220.61120), X64 RyuJIT


|                     Method |  _n |      Mean |     Error |    StdDev |    Median | Ratio | RatioSD |
|--------------------------- |---- |----------:|----------:|----------:|----------:|------:|--------:|
|                   CalcAMix |  50 | 10.873 us | 0.2089 us | 0.1851 us | 10.782 us |  0.82 |    0.01 |
|             CalcAMixArrays |  50 | 13.279 us | 0.0926 us | 0.0723 us | 13.258 us |  1.00 |    0.00 |
|            CalcAMixArrays2 |  50 |  9.492 us | 0.2571 us | 0.7377 us |  9.857 us |  0.70 |    0.05 |
|            CalcAMixArrays3 |  50 |  9.632 us | 0.1921 us | 0.4010 us |  9.436 us |  0.73 |    0.04 |
|              CalcAMixSpans |  50 | 10.204 us | 0.7878 us | 2.1298 us |  9.487 us |  0.81 |    0.20 |
|             CalcAMixSpans2 |  50 |  7.291 us | 0.0284 us | 0.0237 us |  7.296 us |  0.55 |    0.00 |
|             CalcAMixSpans3 |  50 |  5.719 us | 0.1088 us | 0.1252 us |  5.729 us |  0.43 |    0.01 |
|     CalcAMixMirroredArrays |  50 | 12.492 us | 0.2496 us | 0.3959 us | 12.351 us |  0.94 |    0.03 |
|    CalcAMixMirroredArrays2 |  50 |  7.602 us | 0.1500 us | 0.3566 us |  7.728 us |  0.57 |    0.03 |
|    CalcAMixMirroredArrays3 |  50 |  6.095 us | 0.0261 us | 0.0204 us |  6.096 us |  0.46 |    0.00 |
|      CalcAMixMirroredSpans |  50 | 12.993 us | 0.2596 us | 0.4042 us | 13.012 us |  0.97 |    0.03 |
|     CalcAMixMirroredSpans2 |  50 |  7.999 us | 0.2216 us | 0.6065 us |  7.917 us |  0.59 |    0.05 |
|     CalcAMixMirroredSpans3 |  50 |  6.611 us | 0.1343 us | 0.3810 us |  6.682 us |  0.49 |    0.03 |
|    CalcAMixMirroredSpans3a |  50 |  4.233 us | 0.1086 us | 0.3045 us |  4.218 us |  0.30 |    0.01 |
| CalcAMixMirroredPointers3a |  50 |  4.236 us | 0.0846 us | 0.2316 us |  4.173 us |  0.33 |    0.02 |
|                            |     |           |           |           |           |       |         |
|                   CalcAMix | 100 | 48.963 us | 0.9767 us | 1.8345 us | 48.129 us |  0.90 |    0.04 |
|             CalcAMixArrays | 100 | 54.326 us | 1.0678 us | 1.5651 us | 53.187 us |  1.00 |    0.00 |
|            CalcAMixArrays2 | 100 | 35.488 us | 0.7020 us | 1.5112 us | 35.841 us |  0.65 |    0.03 |
|            CalcAMixArrays3 | 100 | 38.967 us | 0.7766 us | 0.8309 us | 38.692 us |  0.72 |    0.02 |
|              CalcAMixSpans | 100 | 35.771 us | 0.7118 us | 1.8999 us | 36.377 us |  0.66 |    0.04 |
|             CalcAMixSpans2 | 100 | 30.509 us | 0.6017 us | 0.9188 us | 29.996 us |  0.56 |    0.02 |
|             CalcAMixSpans3 | 100 | 22.988 us | 0.4522 us | 0.7173 us | 22.457 us |  0.42 |    0.02 |
|     CalcAMixMirroredArrays | 100 | 51.271 us | 1.0005 us | 1.6989 us | 50.254 us |  0.95 |    0.04 |
|    CalcAMixMirroredArrays2 | 100 | 30.479 us | 0.8260 us | 2.1759 us | 30.657 us |  0.55 |    0.03 |
|    CalcAMixMirroredArrays3 | 100 | 26.234 us | 0.5192 us | 0.8530 us | 26.558 us |  0.48 |    0.02 |
|      CalcAMixMirroredSpans | 100 | 52.468 us | 1.0319 us | 1.1469 us | 51.733 us |  0.97 |    0.04 |
|     CalcAMixMirroredSpans2 | 100 | 30.133 us | 0.6320 us | 1.8335 us | 29.517 us |  0.55 |    0.04 |
|     CalcAMixMirroredSpans3 | 100 | 25.254 us | 0.5472 us | 1.6048 us | 24.164 us |  0.47 |    0.04 |
|    CalcAMixMirroredSpans3a | 100 | 15.743 us | 0.3161 us | 0.9120 us | 15.135 us |  0.29 |    0.02 |
| CalcAMixMirroredPointers3a | 100 | 15.988 us | 0.3186 us | 0.7754 us | 15.551 us |  0.30 |    0.02 |

An interesting result here is that mirroring half of the matrix does seem to have some benefits unlike what seemed to be my results from my post #29. I think that with smaller matrices, the entire matrix fits within the CPU cache, and so doing half the work actually pays off, unlike my huge 10000x10000 matrices when I was just benchmarking just filling vs. mirroring.

Spoiler

using System;
using System.Linq;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Running;

namespace AMixBenchmark
{
    public static class ArrayMaker
    {
        static Random _random = new Random(123);

        public static double[] MakeRandomOneDimensionalArray(int n)
            => Enumerable.Range(0, n)
                         .Select(i => _random.NextDouble())
                         .ToArray();

        public static double[,] MakeRandomTwoDimensionalArray(int n)
        {
            var arr = new double[n, n];
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    arr[i, j] = _random.NextDouble();
            return arr;
        }
    }

    public class Bench
    {
        [Params(50, 100)]
        public int _n;
        double [] _xy;
        double [] _ai;
        double [,] _kijm;
        double [,] _aij;
        double [] _sumai;

        [GlobalSetup]
        public void Setup()
        {
            _xy = ArrayMaker.MakeRandomOneDimensionalArray(_n);
            _ai = ArrayMaker.MakeRandomOneDimensionalArray(_n);
            _kijm = ArrayMaker.MakeRandomTwoDimensionalArray(_n);
            _aij = new double[_n, _n];
            _sumai = new double[_n];
        }

        [Benchmark] public double CalcAMix() => CalcAmixInner(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixArrays() => CalcAmixArraysInner(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixArrays2() => CalcAmixArraysInner2(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixArrays3() => CalcAmixArraysInner3(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixSpans() => CalcAmixSpansInner(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixSpans2() => CalcAmixSpansInner2(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixSpans3() => CalcAmixSpansInner3(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredArrays() => CalcAmixMirroredArraysInner(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredArrays2() => CalcAmixMirroredArraysInner2(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredArrays3() => CalcAmixMirroredArraysInner3(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredSpans() => CalcAmixMirroredSpansInner(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredSpans2() => CalcAmixMirroredSpansInner2(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredSpans3() => CalcAmixMirroredSpansInner3(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredSpans3a() => CalcAmixMirroredSpansInner3a(_n, _xy, _ai, _kijm, _aij, _sumai);
        [Benchmark] public double CalcAMixMirroredPointers3a() => CalcAmixMirroredPointersInner3a(_n, _xy, _ai, _kijm, _aij, _sumai);

        static public unsafe double CalcAmixInner(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            double amix = 0;
            int UpperBound0 = Aij.GetUpperBound(0);
            int UpperBound1 = Aij.GetUpperBound(1);

            fixed (double* ptr = Aij, aiPtr = ai, kijmPtr = kijm)
            {
                double* element = ptr;
                double* kijelement = kijmPtr;
                double* aielement = aiPtr;
                for (int i = 0; i <= UpperBound0; i++)
                {
                    for (int j = 0; j <= UpperBound1; j++)
                    {
                        if (i == j)
                            *element = *aielement;
                        else
                            *element = Math.Sqrt(*(aiPtr + i) * *(aiPtr + j)) * (1 - *kijelement);

                        element++;
                        kijelement++;
                    }

                    aielement++;
                }
            }

            fixed (double* ptrAij = Aij, ptrXY = XY, ptrsumAI = sumAi)
            {
                double* elementAij = ptrAij;
                double* elementsumAI = ptrsumAI;
                double* elementXy = ptrXY;

                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        *elementsumAI += *elementXy * *elementAij;
                        elementXy++;
                        elementAij++;
                    }

                    amix += *elementsumAI * *(ptrXY + i);
                    elementXy = ptrXY;
                    elementsumAI++;
                }
            }
            return amix;
        }

        static public double CalcAmixArraysInner(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i == j)
                        Aij[i, j] = ai[i];
                    else
                        Aij[i, j] = Math.Sqrt(ai[i] * ai[j]) * (1 - kijm[i, j]);
                }
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    sumAi[i] += XY[j] * Aij[i, j];

                amix += sumAi[i] * XY[i];
            }
            return amix;
        }

        static public double CalcAmixArraysInner2(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                for (int j = 0; j < n; j++)
                    Aij[i, j] = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);

                Aij[i, i] = aii;
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                for (int j = 0; j < n; j++)
                    sum += XY[j] * Aij[i, j];

                sumAi[i] = sum;
                amix += sum * XY[i];
            }
            return amix;
        }

        static public double CalcAmixArraysInner3(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            double amix = 0;

            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                double sum = XY[i] * aii;

                for (int j = 0; j < n; j++)
                {
                    var value = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                    Aij[i, j] = value;
                    sum += XY[j] * value;
                }

                sum -= XY[i] * Aij[i, i];
                Aij[i, i] = aii;

                sumAi[i] = sum;
                amix += sum * XY[i];
            }

            return amix;
        }

        static public double CalcAmixSpansInner(int n,
                                                ReadOnlySpan<double> XY,
                                                ReadOnlySpan<double> ai,
                                                double[,] kijm,
                                                double[,] Aij,
                                                Span<double> sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                for (int j = 0; j < n; j++)
                    Aij[i, j] = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);

                Aij[i, i] = aii;
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                for (int j = 0; j < n; j++)
                    sum += XY[j] * Aij[i, j];

                sumAi[i] = sum;
                amix += sum * XY[i];
            }
            return amix;
        }

        static public unsafe double CalcAmixSpansInner2(int n,
                                                 ReadOnlySpan<double> XY,
                                                 ReadOnlySpan<double> ai,
                                                 double[,] kijm,
                                                 double[,] Aij,
                                                 Span<double> sumAi)
        {
            fixed (double* kijmPtr = kijm, AijPtr = Aij)
            {
                var kijmSpan = new ReadOnlySpan<double>(kijmPtr, kijm.Length);
                var kijmRowISpan = kijmSpan.Slice(0);
                var AijSpan = new Span<double>(AijPtr, Aij.Length);
                var AijRowISpan = AijSpan.Slice(0);

                for (int i = 0; i < n; i++)
                {
                    var aii = ai[i];
                    for (int j = 0; j < n; j++)
                        AijRowISpan[j] = Math.Sqrt(aii * ai[j]) * (1 - kijmRowISpan[j]);

                    AijRowISpan[i] = aii;
                    kijmRowISpan = kijmRowISpan.Slice(n);
                    AijRowISpan = AijRowISpan.Slice(n);
                }

                double amix = 0;
                AijRowISpan = AijSpan.Slice(0);
                for (int i = 0; i < n; i++)
                {
                    double sum = 0;
                    for (int j = 0; j < n; j++)
                        sum += XY[j] * AijRowISpan[j];

                    sumAi[i] = sum;
                    amix += sum * XY[i];
                    AijRowISpan = AijRowISpan.Slice(n);
                }
                return amix;
            }
        }

        static public unsafe double CalcAmixSpansInner3(int n,
                                                        ReadOnlySpan<double> XY,
                                                        ReadOnlySpan<double> ai,
                                                        double[,] kijm,
                                                        double[,] Aij,
                                                        Span<double> sumAi)
        {
            fixed (double* kijmPtr = kijm, AijPtr = Aij)
            {
                var kijmSpan = new ReadOnlySpan<double>(kijmPtr, kijm.Length);
                var kijmRowISpan = kijmSpan.Slice(0);
                var AijSpan = new Span<double>(AijPtr, Aij.Length);
                var AijRowISpan = AijSpan.Slice(0);
                double amix = 0;

                for (int i = 0; i < n; i++)
                {
                    var aii = ai[i];
                    double sum = XY[i] * aii;

                    for (int j = 0; j < n; j++)
                    {
                        var value = Math.Sqrt(aii * ai[j]) * (1 - kijmRowISpan[j]);
                        AijRowISpan[j] = value;
                        sum += XY[j] * value;
                    }

                    sum -= XY[i] * AijRowISpan[i];
                    AijRowISpan[i] = aii;

                    sumAi[i] = sum;
                    amix += sum * XY[i];

                    kijmRowISpan = kijmRowISpan.Slice(n);
                    AijRowISpan = AijRowISpan.Slice(n);
                }

                return amix;
            }
        }

        static public double CalcAmixMirroredArraysInner(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                Aij[i, i] = ai[i];
                for (int j = i + 1; j < n; j++)
                {
                    Aij[i, j] = Math.Sqrt(ai[i] * ai[j]) * (1 - kijm[i, j]);
                    Aij[j, i] = Aij[i, j];
                }
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    sumAi[i] += XY[j] * Aij[i, j];

                amix += sumAi[i] * XY[i];
            }
            return amix;
        }

        static public double CalcAmixMirroredArraysInner2(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                Aij[i, i] = aii;
                for (int j = i + 1; j < n; j++)
                {
                    var value = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                    Aij[i, j] = value;
                    Aij[j, i] = value;
                }
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                for (int j = 0; j < n; j++)
                    sum += XY[j] * Aij[i, j];

                sumAi[i] = sum;
                amix += sum * XY[i];
            }
            return amix;
        }

        static public double CalcAmixMirroredArraysInner3(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            double amix = 0;

            for (int i = 0; i < n; i++)
            {
                double sum = 0;

                for (int j = 0; j < i; j++)
                    sum += XY[j] * Aij[i, j];

                var aii = ai[i];
                Aij[i, i] = aii;
                sum += XY[i] * aii;

                for (int j = i + 1; j < n; j++)
                {
                    var value = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                    Aij[i, j] = value;
                    Aij[j, i] = value;
                    sum += XY[j] * value;
                }

                sumAi[i] = sum;
                amix += sum * XY[i];
            }

            return amix;
        }

        static public double CalcAmixMirroredSpansInner(int n,
                                                        ReadOnlySpan<double> XY,
                                                        ReadOnlySpan<double> ai,
                                                        double[,] kijm,
                                                        double[,] Aij,
                                                        Span<double> sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                Aij[i, i] = ai[i];
                for (int j = i + 1; j < n; j++)
                {
                    Aij[i, j] = Math.Sqrt(ai[i] * ai[j]) * (1 - kijm[i, j]);
                    Aij[j, i] = Aij[i, j];
                }
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    sumAi[i] += XY[j] * Aij[i, j];

                amix += sumAi[i] * XY[i];
            }
            return amix;
        }

        static public double CalcAmixMirroredSpansInner2(int n,
                                                         ReadOnlySpan<double> XY,
                                                         ReadOnlySpan<double> ai,
                                                         double[,] kijm,
                                                         double[,] Aij,
                                                         Span<double> sumAi)
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                Aij[i, i] = aii;
                for (int j = i + 1; j < n; j++)
                {
                    var value = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                    Aij[i, j] = value;
                    Aij[j, i] = value;
                }
            }

            double amix = 0;
            for (int i = 0; i < n; i++)
            {
                double sum = 0;
                for (int j = 0; j < n; j++)
                    sum += XY[j] * Aij[i, j];

                sumAi[i] = sum;
                amix += sum * XY[i];
            }
            return amix;
        }

        static public double CalcAmixMirroredSpansInner3(int n,
                                                         ReadOnlySpan<double> XY,
                                                         ReadOnlySpan<double> ai,
                                                         double[,] kijm,
                                                         double[,] Aij,
                                                         Span<double> sumAi)
        {
            double amix = 0;

            for (int i = 0; i < n; i++)
            {
                double sum = 0;

                for (int j = 0; j < i; j++)
                    sum += XY[j] * Aij[i, j];

                var aii = ai[i];
                Aij[i, i] = aii;
                sum += XY[i] * aii;

                for (int j = i + 1; j < n; j++)
                {
                    var value = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                    Aij[i, j] = value;
                    Aij[j, i] = value;
                    sum += XY[j] * value;
                }

                sumAi[i] = sum;
                amix += sum * XY[i];
            }

            return amix;
        }

        static public unsafe double CalcAmixMirroredSpansInner3a(int n,
                                                                 ReadOnlySpan<double> XY,
                                                                 ReadOnlySpan<double> ai,
                                                                 double[,] kijm,
                                                                 double[,] Aij,
                                                                 Span<double> sumAi)
        {
            fixed (double* kijmPtr = kijm, AijPtr = Aij)
            {
                var kijmSpan = new ReadOnlySpan<double>(kijmPtr, kijm.Length);
                var kijmRowISpan = kijmSpan.Slice(0);
                var AijSpan = new Span<double>(AijPtr, Aij.Length);
                var AijRowISpan = AijSpan.Slice(0);
                double amix = 0;

                for (int i = 0; i < n; i++)
                {
                    double sum = 0;

                    for (int j = 0; j < i; j++)
                        sum += XY[j] * AijRowISpan[j];

                    var aii = ai[i];
                    AijRowISpan[i] = aii;
                    sum += XY[i] * aii;

                    if (i + 1 < n)
                    {
                        var AijColJSpan = AijSpan.Slice((i + 1) * n);
                        for (int j = i + 1; j < n; j++)
                        {
                            var value = Math.Sqrt(aii * ai[j]) * (1 - kijmRowISpan[j]);
                            AijRowISpan[j] = value;
                            AijColJSpan[i] = value;
                            sum += XY[j] * value;

                            AijColJSpan = AijColJSpan.Slice(n);
                        }
                    }

                    sumAi[i] = sum;
                    amix += sum * XY[i];

                    kijmRowISpan = kijmRowISpan.Slice(n);
                    AijRowISpan = AijRowISpan.Slice(n);
                }

                return amix;
            }
        }

        static public unsafe double CalcAmixMirroredPointersInner3a(int n, double[] XY, double[] ai, double[,] kijm, double[,] Aij, double[] sumAi)
        {
            fixed (double* XYPtr = XY, aiPtr = ai, kijmPtr = kijm, AijPtr = Aij, sumAiPtr = sumAi)
            {
                var kijmRowI = kijmPtr;
                var AijDst = AijPtr;
                var sumAiDst = sumAiPtr;
                double amix = 0;

                for (int i = 0; i < n; i++)
                {
                    var XYSrc = XYPtr;
                    double sum = 0;

                    for (int j = 0; j < i; j++)
                        sum += *XYSrc++ * *AijDst++;

                    var aii = aiPtr[i];
                    *AijDst++ = aii;
                    sum += *XYSrc++ * aii;

                    var AijColJSpan = &AijPtr[(i + 1) * n + i];
                    var kijmSrc = &kijmRowI[i + 1];
                    for (int j = i + 1; j < n; j++)
                    {
                        var value = Math.Sqrt(aii * aiPtr[j]) * (1 - *kijmSrc++);
                        *AijDst++ = value;
                        *AijColJSpan = value;
                        sum += *XYSrc++ * value;

                        AijColJSpan += n;
                    }

                    *sumAiDst++ = sum;
                    amix += sum * XYPtr[i];

                    kijmRowI += n;
                }

                return amix;
            }
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            var summary = BenchmarkRunner.Run<Bench>();
        }
    }
}