Array Acces Speed

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
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

Array Access Speed Test:
 [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());
       
        }
 
Solution
Yes your right, it looks odd now you mention it, when i tested it the results looked ok, late here so ill check it tommorrow. If the results are actually correct maybe that routine can be simplified, the final matrix should be symmetrical around the diagonal, some odd thing in the math maybe. Could just be wrong though. Should be able to cut the sqrt calls in half actually, as its symmetrical, but i was too knackered to get it working today. Need to think about it with a fresh head on.

parralell might be interesting!

i can’t reduce the number of calls to this function for pure tower simulations but definately can for streams and general flash calculations.

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
You still get Spans 3 to be the fastest. I dont get that, different processor maybe? (AMD Ryzen 9 3900X). Anyway I agree the arrays version is simpler and easier to maintain and easy to apply the same methodology to my other arrays. Adding a few pointers makes it less readable.

Amazing what a relatively small change made, program is getting pretty decent speedwise.
 

Skydiver

Staff member
Joined
Apr 6, 2019
Messages
5,085
Location
Chesapeake, VA
Programming Experience
10+
As a follow-up, the only way I could get the "back to arrays 3" to run fastest is to run the Debug build instead of the Release build. This was consistent between my laptop (Ryzen 7), my old PC (AMD Phenom II), and my work laptop (iCore 7). Please always do your timings with Release builds.

Also I saw in your OP that you were running these as unit tests. Don't do that because the VS test engine throws all kinds of overhead. Hopefully you've switched over to using the standalone Console app like I had above.

I should really start using BenchmarkDotNet to guarantee better perf testing because it compile a Release build separately and then time that release build. That's among other good things that BenchmarkDotNet does. The downside is that I can't just copy and paste code anymore for others to try without having them first install the BenchmarkDotNet Nuget package, and actually benchmarking will take longer. Maybe I can do best of both worlds -- use my quick TimeIt() code for quick tests, and then get the full benchmark using BenchmarkDotNet.
 

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
All were running in realease mode, I just tried your original code in a console and still get the Span 3 being slower than Back To Arrays 3. Anyway i think maybe weve optimised this routine as far as it can realtistically go, i have some other areas to re-visit, once again thanks for your valuable help!

2021-01-15.png
 

JohnH

C# Forum Moderator
Staff member
Joined
Apr 23, 2011
Messages
1,415
Location
Norway
Programming Experience
10+
That's odd, these are my results running @Skydiver 's test code from post 16 in Release build on a AMD Ryzen 3600X:
C#:
Original pointers: Average: 503 ms
Back to arrays: Average: 592 ms
Back to arrays 2: Average: 424 ms
Back to arrays 3: Average: 325 ms
Spans: Average: 464 ms
Spans 2: Average: 300 ms
Spans 3: Average: 232 ms
 

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
That's odd, these are my results running @Skydiver 's test code from post 16 in Release build on a AMD Ryzen 3600X:
C#:
Original pointers: Average: 503 ms
Back to arrays: Average: 592 ms
Back to arrays 2: Average: 424 ms
Back to arrays 3: Average: 325 ms
Spans: Average: 464 ms
Spans 2: Average: 300 ms
Spans 3: Average: 232 ms
I think i got it, your running in x64 mode?
 

Skydiver

Staff member
Joined
Apr 6, 2019
Messages
5,085
Location
Chesapeake, VA
Programming Experience
10+
Yes. If you want the most performance or of your 64 bit OS running on a 64 bit machine, you should run as a 64 bit process. If you are running as a 32 but process, the 64 bit OS has to pay overhead costs for thunking to and from a child 32 not process.
 

JohnH

C# Forum Moderator
Staff member
Joined
Apr 23, 2011
Messages
1,415
Location
Norway
Programming Experience
10+
I get similar results for all AnyCPU (both x64 and prefer 32 bit), x64 and x86.
 

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
Yes. If you want the most performance or of your 64 bit OS running on a 64 bit machine, you should run as a 64 bit process. If you are running as a 32 but process, the 64 bit OS has to pay overhead costs for thunking to and from a child 32 not process.
Yeah i had it left as Any Cpu, for no particular reason. It’s all a bit faster now.
 

Skydiver

Staff member
Joined
Apr 6, 2019
Messages
5,085
Location
Chesapeake, VA
Programming Experience
10+
What an interesting result... Trying to do just the top diagonal of the matrix and then trying to mirror it to the lower diagonal seems to be slower than just sequentially computing all the matrix elements. I would have never expected doing half the work to actually take longer. Goes to show how working with the CPU cache instead of against it makes a difference.

C#:
BenchmarkDotNet=v0.12.1, OS=Windows 10.0.19041.685 (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 |     Mean |    Error |   StdDev |   Median |
|------------------------------ |---------:|---------:|---------:|---------:|
| DoAllItemsDiagonalSpecialCase | 532.2 ms | 10.55 ms | 29.77 ms | 540.6 ms |
|           DoAllItemsSegmented | 462.0 ms |  9.22 ms | 13.80 ms | 464.9 ms |
|        DoAllItemsDiagonalLast | 499.7 ms |  9.93 ms | 15.47 ms | 498.9 ms |
|        MirrorItemsConditional | 745.8 ms | 14.39 ms | 21.55 ms | 744.3 ms |
|      MirrorItemsDiagonalFirst | 559.7 ms | 10.97 ms | 10.26 ms | 563.0 ms |
|       MirrorItemsDiagonalLast | 564.0 ms | 11.19 ms | 14.16 ms | 563.1 ms |
|                    DoJustHalf | 247.5 ms |  5.29 ms | 15.43 ms | 246.7 ms |
|  DoJustHalfThenMirrorByColumn | 742.9 ms | 14.71 ms | 17.52 ms | 741.0 ms |
|     DoJustHalfThenMirrorByRow | 778.9 ms | 12.73 ms | 27.66 ms | 773.8 ms |

That was from benchmarking the following code:
C#:
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 BenchArray
    {
        int n = 10000;
        double[] ai;
        double[,] kijm;
        double[,] Aij;

        public BenchArray()
        {
            ai = ArrayMaker.MakeRandomOneDimensionalArray(n);
            kijm = ArrayMaker.MakeRandomTwoDimensionalArray(n);
            Aij = new double[n, n];
        }

        [Benchmark]
        public void DoAllItemsDiagonalSpecialCase()
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                for (int j = 0; j < n; j++)
                {
                    if (i == j)
                        Aij[i, j] = aii;
                    else
                        Aij[i, j] = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                }
            }
        }

        [Benchmark]
        public void DoAllItemsSegmented()
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                for (int j = 0; j < i; j++)
                    Aij[i, j] = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);

                Aij[i, i] = aii;

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

        [Benchmark]
        public void DoAllItemsDiagonalLast()
        {
            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;
            }
        }

        [Benchmark]
        public void MirrorItemsConditional()
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                for (int j = 0; j < n; j++)
                {
                    if (i > j)
                        Aij[j, i] = Aij[i, j];
                    else
                        Aij[i, j] = Math.Sqrt(aii * ai[j]) * (1 - kijm[i, j]);
                }
                Aij[i, i] = aii;
            }
        }

        [Benchmark]
        public void MirrorItemsDiagonalFirst()
        {
            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;
                }
            }
        }

        [Benchmark]
        public void MirrorItemsDiagonalLast()
        {
            for (int i = 0; i < n; i++)
            {
                var aii = ai[i];
                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;
                }
                Aij[i, i] = aii;
            }
        }

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

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

            for (int i = 0; i < n; i++)
            {
                for (int j = i + 1; j < n; j++)
                    Aij[j, i] = Aij[i, j];
            }
        }

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

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < i; j++)
                    Aij[i, j] = Aij[j, i];
            }
        }
    }

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

SilverShaded

Well-known member
Joined
Mar 7, 2020
Messages
65
Programming Experience
10+
So now i'm running x64 (release mode) I found that the size of the arrays makes a difference in the relative speed. For large arrays (10000 x 10000) the modified Back to Arrays 3 method is the fastest and Span 3 is a close second.

However, in my particular case the array sizes will be relatively small (but called many thousands of times), array size is determined by the number of trays in a tower which around 50 would be a kind of median (tallest towers would have maybe 100 trays, beyond that they may be split into two physical towers, these are not so common).

If i run the test with an array size of 50 I find "Back to Arrays 3" & "Modifed back to arrays 3 (or 3a)" are equal in performance and the fastest.

Implication being that the best appproach for this particular problem was optimise the array handling using 'simple' C# code (Skydiver thanks for doing that!) and no need for spans or pointers! There are other arrays were i can try the same re-arranging of the code.

The program is now solving 8 distillation towers and a simple flowsheet in about 3.9s, which is PDG!

Now i can move on to re-write the basic flash calculations elsewhere to be 'inside out' (and hence have less calls to the complex thermodynamics routines).
 

Skydiver

Staff member
Joined
Apr 6, 2019
Messages
5,085
Location
Chesapeake, VA
Programming Experience
10+
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:
C#:
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.

C#:
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>();
        }
    }
}
 
Top Bottom