Matrix API support - start with formats?

[jmh530 via Digitalmars-d]

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);
}