www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - Matrix API support - start with formats?

reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
I stumbled upon https://software.intel.com/en-us/node/471374 which gives 
good detail on Intel's Math Kernel Library's data formats for sparse 
matrices.

No doubt other popular linear algebra libraries have similar 
documentation. I was thinking we could start with adding these layouts 
to std, along with a few simple primitives (construction, element/slice 
access, stride etc). Then, people may just use those as they are or link 
with the linalg libraries for specific computations.


Thoughts?

Andrei
Aug 14 2015
next sibling parent Rikki Cattermole <alphaglosined gmail.com> writes:
On 15/08/2015 2:57 a.m., Andrei Alexandrescu wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 which gives
 good detail on Intel's Math Kernel Library's data formats for sparse
 matrices.

 No doubt other popular linear algebra libraries have similar
 documentation. I was thinking we could start with adding these layouts
 to std, along with a few simple primitives (construction, element/slice
 access, stride etc). Then, people may just use those as they are or link
 with the linalg libraries for specific computations.


 Thoughts?

 Andrei


What would we gain other gl3n in Phobos?
Aug 14 2015
prev sibling next sibling parent reply "ponce" <contact gam3sfrommars.fr> writes:
On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


Are sparse matrices a common scenario?

If anything I think small vectors, non-sparse matrixes and 
rectangles/AABB could be part of Phobos before that.
Aug 14 2015
next sibling parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 8/14/15 11:11 AM, ponce wrote:

 On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 which
 gives good detail on Intel's Math Kernel Library's data formats for
 sparse matrices.

 No doubt other popular linear algebra libraries have similar
 documentation. I was thinking we could start with adding these layouts
 to std, along with a few simple primitives (construction,
 element/slice access, stride etc). Then, people may just use those as
 they are or link with the linalg libraries for specific computations.


 Thoughts?

 Andrei


 Are sparse matrices a common scenario?


In machine learning, yes. But order is not as important as the basic 
strategy: layouts + simple/naive primitives in Phobos, elaborate 
primitives in specialized external libraries. -- Andrei
Aug 14 2015
prev sibling next sibling parent "Tofu Ninja" <emmons0 purdue.edu> writes:
On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
 wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


 Are sparse matrices a common scenario?

 If anything I think small vectors, non-sparse matrixes and 
 rectangles/AABB could be part of Phobos before that.


+1
Aug 14 2015
prev sibling parent reply "David Nadlinger" <code klickverbot.at> writes:
On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 Are sparse matrices a common scenario?


Yes. They tend to pop up in virtually all "serious" numerical 
problems in science and engineering, particularly when partial 
differential equations are involved.


 If anything I think small vectors, non-sparse matrixes and 
 rectangles/AABB could be part of Phobos before that.


If you just had a go at it from the CG/gamedev perspective, you'd 
probably end up with an API that is entirely unsuitable for 
larger problems. This might not be a big issue (just have two 
separate APIs), but we'd need to make it a conscious decision.

  — David
Aug 14 2015
next sibling parent reply "Tofu Ninja" <emmons0 purdue.edu> writes:
On Friday, 14 August 2015 at 18:51:51 UTC, David Nadlinger wrote:

 On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 Are sparse matrices a common scenario?


 Yes. They tend to pop up in virtually all "serious" numerical 
 problems in science and engineering, particularly when partial 
 differential equations are involved.


 If anything I think small vectors, non-sparse matrixes and 
 rectangles/AABB could be part of Phobos before that.


 If you just had a go at it from the CG/gamedev perspective, 
 you'd probably end up with an API that is entirely unsuitable 
 for larger problems. This might not be a big issue (just have 
 two separate APIs), but we'd need to make it a conscious 
 decision.

  — David


Personally I would prefer small fixed sized vectors & matrices 
with game dev focused api be kept separate from a big/sparse 
matrix api.
Aug 14 2015
parent reply Manu via Digitalmars-d <digitalmars-d puremagic.com> writes:
On 15 August 2015 at 05:11, Tofu Ninja via Digitalmars-d
<digitalmars-d puremagic.com> wrote:

 On Friday, 14 August 2015 at 18:51:51 UTC, David Nadlinger wrote:

 On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 Are sparse matrices a common scenario?



 Yes. They tend to pop up in virtually all "serious" numerical problems in
 science and engineering, particularly when partial differential equations
 are involved.


 If anything I think small vectors, non-sparse matrixes and
 rectangles/AABB could be part of Phobos before that.



 If you just had a go at it from the CG/gamedev perspective, you'd probably
 end up with an API that is entirely unsuitable for larger problems. This
 might not be a big issue (just have two separate APIs), but we'd need to
 make it a conscious decision.

  — David



 Personally I would prefer small fixed sized vectors & matrices with game dev
 focused api be kept separate from a big/sparse matrix api.


I think gamedev's would unanimously agree on this.
Aug 20 2015
parent Xavier Bigand <flamaros.xavier gmail.com> writes:
Le 21/08/2015 02:30, Manu via Digitalmars-d a écrit :

 On 15 August 2015 at 05:11, Tofu Ninja via Digitalmars-d
 <digitalmars-d puremagic.com> wrote:

 On Friday, 14 August 2015 at 18:51:51 UTC, David Nadlinger wrote:

 On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 Are sparse matrices a common scenario?



 Yes. They tend to pop up in virtually all "serious" numerical problems in
 science and engineering, particularly when partial differential equations
 are involved.


 If anything I think small vectors, non-sparse matrixes and
 rectangles/AABB could be part of Phobos before that.



 If you just had a go at it from the CG/gamedev perspective, you'd probably
 end up with an API that is entirely unsuitable for larger problems. This
 might not be a big issue (just have two separate APIs), but we'd need to
 make it a conscious decision.

   — David



 Personally I would prefer small fixed sized vectors & matrices with game dev
 focused api be kept separate from a big/sparse matrix api.


 I think gamedev's would unanimously agree on this.


Those kind of maths aren't used only for games, their is also editing 
software,...

Their is also spacial relations that can be interesting to have in a 
geometry module. Boost do it and respect those rules : 
http://edndoc.esri.com/arcsde/9.0/general_topics/understand_spatial_relations.htm
For my job it help us to solve a lot of precision issues.

Having this in the standard library is important IMO cause sadly their 
is often a lot of bugs/issues when reimplemented handly. A lot of 
articles on Internet are just wrong.
Aug 21 2015
prev sibling next sibling parent "jmh530" <john.michael.hall gmail.com> writes:
On Friday, 14 August 2015 at 18:51:51 UTC, David Nadlinger wrote:

 If anything I think small vectors, non-sparse matrixes and 
 rectangles/AABB could be part of Phobos before that.


 If you just had a go at it from the CG/gamedev perspective, 
 you'd probably end up with an API that is entirely unsuitable 
 for larger problems. This might not be a big issue (just have 
 two separate APIs), but we'd need to make it a conscious 
 decision.


While I'm interested in working with big (500,000+) dense (and 
sometimes smaller sparse) rectangular matrices, my sense of the D 
forum is that there are more CG/gamedev people using D.

Nevertheless, I think it is important to get the general 
structure right in a way that can be specialized to work with 
small matrices, rectangular matrices, strided matrices, or sparse 
matrices. I don't see an issue with prioritizing what users want 
in terms of fleshing out the APIs, so long as the different ways 
people use matrix libraries is in the back of the designers minds 
when they design it.

I was looking over some of the OpenGL matrix libraries recently. 
They have types like mat4 or vec3. I couldn't help but think that 
template constraints (at least for static arrays) would probably 
make it much easier. For instance, if you have a function that 
takes a vec, you can specialize it to do different things 
depending on if the array is small or large (small, unroll loops, 
large use parallelism, maybe with the decision on when stop 
unrolling loops or starting to use parallelism determined by 
AEOS, as in ATLAS).
Aug 14 2015
prev sibling parent "ponce" <contact gam3sfrommars.fr> writes:
On Friday, 14 August 2015 at 18:51:51 UTC, David Nadlinger wrote:

 On Friday, 14 August 2015 at 15:11:39 UTC, ponce wrote:

 Are sparse matrices a common scenario?


 Yes. They tend to pop up in virtually all "serious" numerical 
 problems in science and engineering, particularly when partial 
 differential equations are involved.


 If anything I think small vectors, non-sparse matrixes and 
 rectangles/AABB could be part of Phobos before that.


 If you just had a go at it from the CG/gamedev perspective, 
 you'd probably end up with an API that is entirely unsuitable 
 for larger problems. This might not be a big issue (just have 
 two separate APIs), but we'd need to make it a conscious 
 decision.


One single API would be great in that case. If someone can pull 
it off.
Aug 14 2015
prev sibling next sibling parent reply "bachmeier" <no spam.com> writes:
On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


One concern I have is the choice of MKL, which due to cost and 
license reasons, many developers will not have on all (or even 
any) of their machines.

I don't work with sparse matrices often so I do not know which 
libraries are most popular, but at a minimum I think it is 
necessary to say you can use it with a popular open source 
library. Given the importance of MKL, it would be a bad idea to 
not offer a compatible format, but it would be equally bad to 
focus on only MKL.
Aug 14 2015
parent reply "jmh530" <john.michael.hall gmail.com> writes:
On Friday, 14 August 2015 at 19:19:24 UTC, bachmeier wrote:

 One concern I have is the choice of MKL, which due to cost and 
 license reasons, many developers will not have on all (or even 
 any) of their machines.

 I don't work with sparse matrices often so I do not know which 
 libraries are most popular, but at a minimum I think it is 
 necessary to say you can use it with a popular open source 
 library. Given the importance of MKL, it would be a bad idea to 
 not offer a compatible format, but it would be equally bad to 
 focus on only MKL.


I agree. MKL is expensive. OpenBLAS is supposed have comparable 
performance, more or less, and it is free. Alternately, ATLAS can 
be used to build BLAS on many different systems.

I would also distinguish between the low level API like BLAS and 
OpenBLAS and higher level APIs. Good higher level APIs allow 
drop-in replacement of lower level APIs (e.g. Armadillo).
Aug 14 2015
parent reply "jmh530" <john.michael.hall gmail.com> writes:
On Friday, 14 August 2015 at 20:23:00 UTC, jmh530 wrote:

 I agree. MKL is expensive. OpenBLAS is supposed have comparable 
 performance, more or less, and it is free. Alternately, ATLAS 
 can be used to build BLAS on many different systems.

 I would also distinguish between the low level API like BLAS 
 and OpenBLAS and higher level APIs. Good higher level APIs 
 allow drop-in replacement of lower level APIs (e.g. Armadillo).


It looks like AMD has open-sourced some of their OpenCL math 
libraries as well.

https://github.com/clMathLibraries/clBLAS
https://github.com/clMathLibraries/clRNG
https://github.com/clMathLibraries/clSPARSE
https://github.com/clMathLibraries/clFFT

They also made some available through
https://github.com/flame
Aug 20 2015
parent "jmh530" <john.michael.hall gmail.com> writes:
On Thursday, 20 August 2015 at 21:28:10 UTC, jmh530 wrote:

 https://github.com/clMathLibraries/clBLAS
 https://github.com/clMathLibraries/clRNG
 https://github.com/clMathLibraries/clSPARSE
 https://github.com/clMathLibraries/clFFT

 They also made some available through
 https://github.com/flame


I probably should have also highlighted the BLIS project in the 
FLAME repo.
https://code.google.com/p/blis/
Aug 20 2015
prev sibling next sibling parent "jmh530" <john.michael.hall gmail.com> writes:
On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


I am very excited that there is a greater focus on this.

I have been following this for a little while:
https://github.com/D-Programming-Language/phobos/pull/3397
It already does or intends to cover many things that I think are 
needed: easily looping through every element of a nd-slice, 
looping by dimension (by row or by column), and deep map. I 
wouldn't be surprised if there are language improvements that can 
facilitate this work.

In terms of formatting, I was playing around with a hybrid array 
struct (below) that is basically a static array, but also sort of 
range-ified. I was using the concept of domains from Chapel, 
albeit at an elementary level that could be significantly 
improved.

For that matter, may I suggest that you look over the Chapel 
specification:
http://chapel.cray.com/spec/spec-0.91.pdf
most importantly the domains and arrays sections. From writing 
that hybrid_array code, I think D lends itself well to some of 
these ideas.

Nevertheless, I agree that interoperability with C-based linear 
algebra libraries is important. If anything I said would make it 
more difficult, then ignore it.

import std.traits;
import std.range;
import std.array : array;
import std.stdio : writeln;
import std.algorithm : map;

struct hybrid_array(T : U[N], U, size_t N)
	if ( isStaticArray!T )
{
	T static_array;
	private auto _length = static_array.length;
	auto domain = iota(0, static_array.length);
	
	this(T x)
	{
		static_array = x;
	}
	
	 property bool empty()
     {
         return domain.empty;
     }

      property size_t length()
     {
         return _length;
     }

      property U front()
     {
		assert(!empty);
         return static_array[domain.front];
     }

     void popFront()
     {
		assert(!empty);
         domain.popFront;
		_length--;
     }
	
	 property hybrid_array save()
     {
         return this;
     }
	
	 property U back()
     {
		assert(!empty);
         return static_array[domain.back];
     }

     void popBack()
     {
		assert(!empty);
         domain.popBack;
		_length++;
     }
	
	U opIndex(size_t val)
     {
		assert(!empty);
         return static_array[domain[val]];
     }
}

void main()
{
	int[5] x = [0, 10, 2, 6, 1];
	//auto squares = map!(a => a * a)(x); //doesn't work
	
	auto range = iota(0, x.length).array;
	
	alias h_array = hybrid_array!(int[5]);
	auto y = h_array(x);
	
	writeln( isRandomAccessRange!(typeof(y)) );
	
	auto squares = map!(a => a * a)(y); //success, does work
	writeln(squares);
}
Aug 14 2015
prev sibling next sibling parent Shachar Shemesh <shachar weka.io> writes:
On 14/08/15 17:57, Andrei Alexandrescu wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 which gives
 good detail on Intel's Math Kernel Library's data formats for sparse
 matrices.

 No doubt other popular linear algebra libraries have similar
 documentation. I was thinking we could start with adding these layouts
 to std, along with a few simple primitives (construction, element/slice
 access, stride etc). Then, people may just use those as they are or link
 with the linalg libraries for specific computations.


 Thoughts?


Please please please make them templated on the type, with clear 
instructions what the type needs to support in order to work.

I've recently tried to find an implementation of matrix inversion over a 
Z2 field, and found that that NONE of the standard implementations 
supports that.

Shachar
Aug 15 2015
prev sibling next sibling parent "Daniel Shapero" <shapero.daniel gmail.com> writes:
On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
wrote:

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


For comparison, some other sparse linear algebra libraries worth 
looking at are:

* Eigen (http://eigen.tuxfamily.org/)
* PETSc (http://www.mcs.anl.gov/petsc/)
* the Epetra sub-package of Trilinos (https://trilinos.org/)
* OSKI, which aims to do for sparse matrices what ATLAS does for 
dense (http://bebop.cs.berkeley.edu/oski/)

MKL is missing some important sparse matrix formats -- the 
ellpack and jagged diagonal formats, which are very well suited 
for SIMD processors, and the MSR/MSC format, which makes 
multigrid smoothing faster.

I don't know enough to weigh in on whether sparse matrix algebra 
is best left to the individual libraries or put in the language's 
standard library. It's common practice whenever one needs sparse 
matrices to write wrapper classes that can use either PETSc or 
Trilinos, depending on what the users have installed on their 
systems. For example, this is the approach taken in the finite 
element library deal.II (http://dealii.org/). Moreover, the big 
sparse matrix libraries can all use each other; PETSc can link to 
the solvers in Trilinos, both PETSc and Trilinos can link to 
OSKI, etc. If there were some support in the standard library for 
this functionality, it could obviate the need for some of the 
glue code in C/C++.
Aug 15 2015
prev sibling next sibling parent reply "dominik" <dschoerk gmx.at> writes:
coming from a computer graphics and image processing background 
i'd love to have a standard way to deal with matrices.

But matrices come in a lot of different flavors, which need 
different implementations.
- computer graphics applications need mostly 1d-4d vectors and 
matrices.
   encapsulation of transformation types: 
rigid/similarity/affine/linear/projective give the
   opportunity to choose simpler algorithms for e.g. inversion of 
a matrix
- images are usually represented as dense 2d matrix (or 3d if 
multiple colors channels are present)
- equation systems use dense and sparse matrices, and use lots of 
different data structures
   to be optimal for different algorithms. yale dok lil coo to 
name some sparse formats.
   mathematically they have different properties like 
square/rectangular positive/negative
   definit diagonal
- matlab style matrices
   for convenience it should be possible to resize, and 
concatenate matrices in different
   dimensions, e.g. to compose images or to augment equation 
systems

on top of that possibilities to solve those equation systems 
would be handy:
with c++ i mostly use eigen, opennl or ceres

this was just a list of things i thought should be kept in mind 
when implementing a matrix api - the matrix api doesnt have to 
implement all this

imho the matrix api could provide some basic types and a general 
interface to access such matrices. for start a simple dense and 
sparse n-dimensional matrix format with the possibility to 
construct and change matrices. slicing in n-dimensions as simple 
as in matlab would be awesome.

such an api would also be a great base for writing graph based 
algorithms like markov random fields ...

hopes and wishes :)
best dominik
Aug 21 2015
parent "dominik" <dschoerk gmx.at> writes:
On Friday, 21 August 2015 at 11:31:16 UTC, dominik wrote:

 - computer graphics applications need mostly 1d-4d vectors and 
 matrices.


vec2, vec3, vec4, mat2, mat3, mat4 was what i ment obviously
Aug 21 2015
prev sibling parent reply "ZombineDev" <valid_email he.re> writes:
On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


One thing that will make D really shine is to implement something 
like this:
https://www.youtube.com/watch?v=hfn0BVOegac

Since Blaze [1] is open source all we need to do is to provide D 
wrappers on their highly optimized kernels. My intuition is that 
by using D's generative features we may be able implement this 
with significantly less effort than with C++. Their repository 
even contains benchmarks which we can use to verify that our 
wrappers won't incur overhead compared to C++.

[1]: https://bitbucket.org/blaze-lib/blaze
Sep 01 2015
next sibling parent "Laeeth Isharc" <laeethnospam nospamlaeeth.com> writes:
On Tuesday, 1 September 2015 at 14:00:52 UTC, ZombineDev wrote:

 On Friday, 14 August 2015 at 14:57:19 UTC, Andrei Alexandrescu 
 wrote:

 I stumbled upon https://software.intel.com/en-us/node/471374 
 which gives good detail on Intel's Math Kernel Library's data 
 formats for sparse matrices.

 No doubt other popular linear algebra libraries have similar 
 documentation. I was thinking we could start with adding these 
 layouts to std, along with a few simple primitives 
 (construction, element/slice access, stride etc). Then, people 
 may just use those as they are or link with the linalg 
 libraries for specific computations.


 Thoughts?

 Andrei


 One thing that will make D really shine is to implement 
 something like this:
 https://www.youtube.com/watch?v=hfn0BVOegac

 Since Blaze [1] is open source all we need to do is to provide 
 D wrappers on their highly optimized kernels. My intuition is 
 that by using D's generative features we may be able implement 
 this with significantly less effort than with C++. Their 
 repository even contains benchmarks which we can use to verify 
 that our wrappers won't incur overhead compared to C++.

 [1]: https://bitbucket.org/blaze-lib/blaze


I was looking at blaze the other day and wondering just the same. 
  I haven't used it, and don't have a great feeling for what would 
be involved.  Maybe John Colvin does.
Sep 02 2015
prev sibling parent "jmh530" <john.michael.hall gmail.com> writes:
On Tuesday, 1 September 2015 at 14:00:52 UTC, ZombineDev wrote:

 One thing that will make D really shine is to implement 
 something like this:
 https://www.youtube.com/watch?v=hfn0BVOegac



Thanks for posting that. It provides a nuanced discussion of 
expression templates. The final conclusion is that C++ should be 
very careful in choosing a linear algebra library for the 
standard. I think that applies to D just as well.
Sep 02 2015