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

Reply via email to