Hi All,

The Dlang multidimensional range type, ndslice, is a struct composed a an iterator, lengths and possibly strides. It does not own memory and does not know anything about its content. ndslice is a faster and extended version of numpy.ndarray.

After some work on commercial projects based on Lubeck[1] and ndslice I figure out what API and memory management is required to make Dlang super fast and  math friendly in the same time.

The concept is the following:
1. All memory is managed by a global BetterC thread safe ARC allocator. Optionally the allocator can be overloaded.
2. User can take an internal ndslice to use mir.ndslice API internally in functions.
2. auto matrixB = matrixA; // increase ARC
3. auto matrixB = matrixA.dup; // allocates new matrix
4. matrix[i] returns a Vec and increase ARC, matrix[i, j] returns a content of the cell.
5. Clever `=` expression based syntax. For example:

   // performs CBLAS call of GEMM and does zero memory allocations
   C = alpha * A * B + beta * C;

`Mat` and other types will support any numeric types, PODlike structs, plus special overload for `bool` based on `bitwise` [2].

I have a lot of work for next months, but looking for a good opportunity to make Mat happen.

For contributing or co-financing:
Ilya Yaroshenko at
gmail com

Best Regards,
Ilya

[1] https://github.com/kaleidicassociates/lubeck
[2] http://docs.algorithm.dlang.io/latest/mir_ndslice_topology.html#bitwise
[3] http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html#gaeda3cbd99c8fb834a60a6412878226e1

On Tuesday, 13 March 2018 at 03:37:36 UTC, 9il wrote:
> [snip]
> 4. matrix[i] returns a Vec and increase ARC, matrix[i, j] returns a content of the cell.
> 

I'm not 100% sold on matrix[i] returning a Vec instead of a 1-dimensional matrix. R does something similar and you have to convert things back to a matrix for some computations more often than I'd like. If functions can easily take both Mat and Vec types in a relatively painless fashion, then I wouldn't have an issue with it.

> 5. Clever `=` expression based syntax. For example:
>
>    // performs CBLAS call of GEMM and does zero memory allocations
>    C = alpha * A * B + beta * C;

You might want to explain this in more detail. I saw expression and my head went to expression templates, but that doesn't seem to be what you're talking about (overloading opAssign?)

> I have a lot of work for next months, but looking for a good opportunity to make Mat happen.
>

+1

With respect to the title, the benefit of special matrix types is when we can call functions (lapack or otherwise) that are optimized for those types. If you want the best performance for mir, then I think that's what it would take. I'm not sure how much you've thought about this. For instance, I understand from graphics libraries that if you're only working with a particular size matrix (say 3x3), then you can generate faster code than if you're working with general matrices.

In addition, performance is not the only thing a new user to mir would care about They likely would also care about ease-of-use [1] and documentation. Hopefully these continue to improve.

What's TMMat?

Diag seems like it would be a special case of sparse matrices, though diag is probably simpler to implement.

[1] Would it be seamless to add a Mat to a Diag? Also what happens to the api when you add 10 different matrix types and need to think about all the interactions.

On Tuesday, 13 March 2018 at 03:37:36 UTC, 9il wrote:
> Hi All,
>
> The Dlang multidimensional range type, ndslice, is a struct composed a an iterator, lengths and possibly strides. It does not own memory and does not know anything about its content. ndslice is a faster and extended version of numpy.ndarray.
>
> After some work on commercial projects based on Lubeck[1] and ndslice I figure out what API and memory management is required to make Dlang super fast and  math friendly in the same time.
>
> The concept is the following:
> 1. All memory is managed by a global BetterC thread safe ARC allocator. Optionally the allocator can be overloaded.
> 2. User can take an internal ndslice to use mir.ndslice API internally in functions.
> 2. auto matrixB = matrixA; // increase ARC
> 3. auto matrixB = matrixA.dup; // allocates new matrix
> 4. matrix[i] returns a Vec and increase ARC, matrix[i, j] returns a content of the cell.
> 5. Clever `=` expression based syntax. For example:
>
>    // performs CBLAS call of GEMM and does zero memory allocations
>    C = alpha * A * B + beta * C;
>
> `Mat` and other types will support any numeric types, PODlike structs, plus special overload for `bool` based on `bitwise` [2].
>
> I have a lot of work for next months, but looking for a good opportunity to make Mat happen.
>
> For contributing or co-financing:
> Ilya Yaroshenko at
> gmail com
>
> Best Regards,
> Ilya
>
> [1] https://github.com/kaleidicassociates/lubeck
> [2] http://docs.algorithm.dlang.io/latest/mir_ndslice_topology.html#bitwise
> [3] http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html#gaeda3cbd99c8fb834a60a6412878226e1

Well if D had:
- a good matrix type, supporting all numeric types (absolutely crucial: including complex!)
- *very important*: with operations syntax corresponding to the well-established standard mathematical conventions (including in labelling the elements!)
- with superfast LU decomposition, det and similar (perhaps via LAPACK)
- able to recognize/carry a tag for (anti)symmetric, hermitian, ... and exploit these in computations
then I'd be more seriously considering switching to D.

Your suggestion [4] that matrix[i] returns a Vec is perhaps too inflexible. What one needs sometimes is to return a row, or a column of a matrix, so a notation like matrix[i, ..] or matrix[.., j] returning respectively a row or column would be useful.

I'd be happy to help out/give advice from the viewpoint of a scientist trying to "graduate" away from C++ without having to sacrifice performance.

On Tuesday, 13 March 2018 at 04:35:53 UTC, jmh530 wrote:
>> 5. Clever `=` expression based syntax. For example:
>>
>>    // performs CBLAS call of GEMM and does zero memory allocations
>>    C = alpha * A * B + beta * C;
>
> You might want to explain this in more detail. I saw expression and my head went to expression templates, but that doesn't seem to be what you're talking about (overloading opAssign?)

Expression templates plus overloading opAssign.

>> I have a lot of work for next months, but looking for a good opportunity to make Mat happen.
>>
>
> +1
>
> With respect to the title, the benefit of special matrix types is when we can call functions (lapack or otherwise) that are optimized for those types. If you want the best performance for mir, then I think that's what it would take. I'm not sure how much you've thought about this. For instance, I understand from graphics libraries that if you're only working with a particular size matrix (say 3x3), then you can generate faster code than if you're working with general matrices.
>
> In addition, performance is not the only thing a new user to mir would care about They likely would also care about ease-of-use [1] and documentation. Hopefully these continue to improve.
>
> What's TMMat?

TMat is a transposed matrix. Not sure for now if it would be required.

>
> Diag seems like it would be a special case of sparse matrices, though diag is probably simpler to implement.

Diag should be an independent type. But not in the first release.

> [1] Would it be seamless to add a Mat to a Diag?

Mat C = M + D; // D - diagonal.
D += M; // Fails at compile time, but
D += M.diag; // pass

> Also what happens to the api when you add 10 different matrix types and need to think about all the interactions.

User API would be simple. But it requires clever compile and runtime logic to do the best with expression templates.

On Tuesday, 13 March 2018 at 05:36:06 UTC, J-S Caux wrote:
>
> Your suggestion [4] that matrix[i] returns a Vec is perhaps too inflexible. What one needs sometimes is to return a row, or a column of a matrix, so a notation like matrix[i, ..] or matrix[.., j] returning respectively a row or column would be useful.

auto row = matrix[i]; // or matrix[i, 0 .. $];
auto col = matrix[0 .. $, j];

On Tuesday, 13 March 2018 at 10:39:29 UTC, 9il wrote:
> On Tuesday, 13 March 2018 at 05:36:06 UTC, J-S Caux wrote:
>>
>> Your suggestion [4] that matrix[i] returns a Vec is perhaps too inflexible. What one needs sometimes is to return a row, or a column of a matrix, so a notation like matrix[i, ..] or matrix[.., j] returning respectively a row or column would be useful.
>
> auto row = matrix[i]; // or matrix[i, 0 .. $];
> auto col = matrix[0 .. $, j];

Some kind of improvement that replaces 0 .. $ with some shorter syntax has been brought up in the past.
https://github.com/libmir/mir-algorithm/issues/53

On Tuesday, 13 March 2018 at 03:37:36 UTC, 9il wrote:
[...]
> 5. Clever `=` expression based syntax. For example:
>
>    // performs CBLAS call of GEMM and does zero memory allocations
>    C = alpha * A * B + beta * C;
>
[...]
My answer is: Yes.

If D with Lubeck would have such a convenient way to write matrix algebra,
using well designed operator overloading, this might attract more users.

Unfortunately my own time using this kind of math daily is long ago,
but I would like to help as a tester.

I will look in my master thesis and see if I may use Lubeck for the calculation done.
I was doing eigenvalue with FEM and in an other project partial differential equation with a matrix based approximation schema... so this may bring my "gray cells" to work again :-)

I was using a C++ lib with operator overloading resulting in quite convenient expressions.
The point that Java has no operator overloading, was the trigger for me to dislike the language.  :-)

On Tuesday, 13 March 2018 at 10:39:29 UTC, 9il wrote:
> On Tuesday, 13 March 2018 at 05:36:06 UTC, J-S Caux wrote:
>>
>> Your suggestion [4] that matrix[i] returns a Vec is perhaps too inflexible. What one needs sometimes is to return a row, or a column of a matrix, so a notation like matrix[i, ..] or matrix[.., j] returning respectively a row or column would be useful.
>
> auto row = matrix[i]; // or matrix[i, 0 .. $];
> auto col = matrix[0 .. $, j];

Some kind of improvement that replaces 0 .. $ with some shorter syntax has been brought up in the past.
https://github.com/libmir/mir-algorithm/issues/53

On Tuesday, 13 March 2018 at 12:16:27 UTC, jmh530 wrote:
>
> Some kind of improvement that replaces 0 .. $ with some shorter syntax has been brought up in the past.
> https://github.com/libmir/mir-algorithm/issues/53

Sorry for double post.

On Tuesday, 13 March 2018 at 10:35:15 UTC, 9il wrote:
> On Tuesday, 13 March 2018 at 04:35:53 UTC, jmh530 wrote:
>> [snip]
>>
>> What's TMMat?
>
> TMat is a transposed matrix. Not sure for now if it would be required.
>

There are some people who like being able to specify a whether a matrix has column or row layout. Would an option to control this be the same thing?