On Thursday, 27 February 2020 at 23:15:28 UTC, p.shkadzko wrote:
> And it works effortlessly!
> Sum of two 5000 x 6000 int arrays is just 0.105 sec! (on a Windows machine though but with weaker CPU).
>
> I bet using mir.ndslice instead of D arrays would be even faster.

Yes, the output for the following benchmark shows that Mir is 43% faster.
However, when I have checked the assembler output, both Mir and Std (really LDC in both cases) generate almost the same and best possible loops with AVX instructions for summation.

In another hand, Mir is faster because it generates random matrixes faster and uses uninitialized memory for the summation target.

Output:
```
std: 426 ms, 432 μs, and 1 hnsec |10
mir: 297 ms, 694 μs, and 3 hnsecs |10
```

Run command:

`dub --build=release --single --compiler=ldc2 test.d`

Note that -mcpu=native flag is passed to LDC.

Source:
```
/+dub.sdl:
dependency "mir-algorithm" version="~>3.7.17"
dependency "mir-random" version="~>2.2.10"
dflags "-mcpu=native" platform="ldc"
+/

int val;

void testStd()
{
    pragma(inline, false);
    static struct Matrix(T)
    {
        import std.range;
        T[] elems;
        int cols;

        T[][] to2D()
        {
            return elems.chunks(cols).array;
        }
    }

    static auto matrixSum(Matrix!int m1, Matrix!int m2)
    {
        Matrix!int m3;
        m3.cols = m1.cols;
        m3.elems.length = m1.elems.length;
        m3.elems[] = m1.elems[] + m2.elems[];
        return m3.to2D;
    }

    static T[] rndArr(T)(in T max, in int elems)
    {
        import std.random;
        import std.range;
        Xorshift rnd;
        return generate(() => uniform(0, max, rnd)).take(elems).array;
    }
    auto m1 = Matrix!int(rndArr!int(10, 5000 * 6000), 6000);
    auto m2 = Matrix!int(rndArr!int(10, 5000 * 6000), 6000);
    auto m3 = matrixSum(m1, m2);
    val = m3[$-1][$-1];
}

void testMir()
{
    pragma(inline, false);
    import mir.ndslice;
    import mir.random: threadLocal;
    import mir.random.variable: uniformVar;
    import mir.random.algorithm: randomSlice;
    import mir.random.engine.xorshift;

    auto m1 = threadLocal!Xorshift.randomSlice(uniformVar!int(0, 10), [5000, 6000]);
    auto m2 = threadLocal!Xorshift.randomSlice(uniformVar!int(0, 10), [5000, 6000]);
    auto m3 = slice(m1 + m2);
    val = m3[$-1][$-1];
}

void main()
{
	import std.datetime.stopwatch;
    import std.stdio;
    import core.memory;
    GC.disable;
    StopWatch clock;
    clock.reset;
    clock.start;
    testStd;
    clock.stop;
    writeln("std: ", clock.peek, " |", val);
    clock.reset;
    clock.start;
    testMir;
    clock.stop;
    writeln("mir: ", clock.peek, " |", val);
}
```

On Thursday, 27 February 2020 at 14:15:26 UTC, p.shkadzko wrote:
> I'd like to sum 2D arrays. Let's create 2 random 2D arrays and sum them.
>
> ```
> import std.random : Xorshift, unpredictableSeed, uniform;
> import std.range : generate, take, chunks;
> import std.array : array;
>
> static T[][] rndMatrix(T)(T max, in int rows, in int cols)
> {
>     Xorshift rnd;
>     rnd.seed(unpredictableSeed);
>     const amount = rows * cols;
>     return generate(() => uniform(0, max, rnd)).take(amount).array.chunks(cols).array;
> }
>
> void main() {
>     int[][] m1 = rndMatrix(10, 2, 3);
>     int[][] m2 = rndMatrix(10, 2, 3);
>
>     auto c = m1[] + m2[];
> }
> ```
>

Maybe this is already clear, but it is important to highlight that 2D arrays and arrays of arrays are two different things.

int[][] is an array of arrays, for each outer index the element is an array that has its own allocated memory and length.

2D arrays are not provided by the language, but can be implemented by defining a type with the required operator overloads.

Using an array of arrays for a 10000x10000 matrix requires 10001 allocations while a dedicated 2D array implementation needs only 1;


Example of an array of arrays where the inner arrays have different lengths:

------------
module test;
void main()
{
    import std.stdio;

    int[][] a;
    a.length=3;

    a[0]=[1,2,3];
    a[1]=[3,4];
    a[2]=[];
    writeln(a);
}
------------


Your Example with a minimal 2D array.


------------
module test2;

import std.random : Xorshift, unpredictableSeed, uniform;
import std.range : generate, take, chunks;
import std.array : array;
import std.stdio : writeln;

struct Matrix(T)
{
 int rows;
 T[] data;
 alias data this;
 int cols() {return cast(int) data.length/rows;}
 this(int r, int c) { data=new int[r*c]; rows=r;}
 this(int r, int c, T[] d) {assert(r*c==data.length); data=d; rows=r; }

 auto opIndex(int r, int c) {return data[rows*c+r];}

}

auto rndMatrix(T)(T max, in int rows, in int cols)
{
    Xorshift rnd;
    rnd.seed(unpredictableSeed);
    const amount = rows * cols;
    return Matrix!T(rows,cols,generate(() => uniform(0, max, rnd)).take(amount).array);
}

void main() {
    auto m1 = rndMatrix(10, 2, 3);
    auto m2 = rndMatrix(10, 2, 3);

    auto c = Matrix!int(2,3);
    c[] = m1[] + m2[];
    writeln(m1[1,2]);
    writeln(m2[1,2]);
    writeln(c[1,2]);
}
----------------

See https://dlang.org/spec/operatoroverloading.html#array-ops for a better overview of the required operators or mir.ndslice for an nD implementation.

On Friday, 28 February 2020 at 16:51:10 UTC, AB wrote:

> See https://dlang.org/spec/operatoroverloading.html#array-ops for a better overview of the required operators or mir.ndslice for an nD implementation.

Here's an old version of some of the things I've been using:

https://bitbucket.org/bachmeil/dmdgretl/src/67a6c5dbf95f23fa753bfd590bc464147cbdc5cc/inst/gretl/matrix.d#lines-307

It has multidimensional slicing and such.

On Friday, 28 February 2020 at 16:51:10 UTC, AB wrote:
> On Thursday, 27 February 2020 at 14:15:26 UTC, p.shkadzko wrote:
>>[...]
>
> Your Example with a minimal 2D array.
>
> ------------
> module test2;
>
> import std.random : Xorshift, unpredictableSeed, uniform;
> import std.range : generate, take, chunks;
> import std.array : array;
> import std.stdio : writeln;
>
> struct Matrix(T)
> {
>  int rows;
>  T[] data;
>  alias data this;
>  int cols() {return cast(int) data.length/rows;}
>  this(int r, int c) { data=new int[r*c]; rows=r;}
>  this(int r, int c, T[] d) {assert(r*c==data.length); data=d; rows=r; }
>
>  auto opIndex(int r, int c) {return data[rows*c+r];}
>
> }

Can you please explain what is the purpose of "alias data this" in your Matrix struct? As I remember "alias <member> this" is used for implicit type conversions but I don't see where "data" is converted.

On Saturday, 29 February 2020 at 19:04:12 UTC, p.shkadzko wrote:
> On Friday, 28 February 2020 at 16:51:10 UTC, AB wrote:
>> On Thursday, 27 February 2020 at 14:15:26 UTC, p.shkadzko wrote:
>>>[...]
>>
>> Your Example with a minimal 2D array.
>>
>> ------------
>> module test2;
>>
>> import std.random : Xorshift, unpredictableSeed, uniform;
>> import std.range : generate, take, chunks;
>> import std.array : array;
>> import std.stdio : writeln;
>>
>> struct Matrix(T)
>> {
>>  int rows;
>>  T[] data;
>>  alias data this;
>>  int cols() {return cast(int) data.length/rows;}
>>  this(int r, int c) { data=new int[r*c]; rows=r;}
>>  this(int r, int c, T[] d) {assert(r*c==data.length); data=d; rows=r; }
>>
>>  auto opIndex(int r, int c) {return data[rows*c+r];}
>>
>> }
>
> Can you please explain what is the purpose of "alias data this" in your Matrix struct? As I remember "alias <member> this" is used for implicit type conversions but I don't see where "data" is converted.

Without "alias data this" the call

    c[] = m1[] + m2[];

should be written as

    c.data[] = m1.data[] + m2.data[];

With "alias data this" if some members of Matrix are not defined, but they are available for Matrix.data, the members of Matrix.data will be used.

It is more about names lookup rules than type conversions.