Linear algebra library

[H. S. Teoh]

On Tue, Apr 09, 2013 at 06:57:08PM +0200, Maksim Zholudev wrote:
> On Tuesday, 9 April 2013 at 14:41:09 UTC, Denis Shelomovskij wrote:
[...]
> >Also, if we are talking about matrices, very sad a solution for
> >enhancement 6798 [1] isn't merged yet.
> >
> >[1] http://d.puremagic.com/issues/show_bug.cgi?id=6798
> 
> Yes, I know about this. At the beginning I tried to implement
> multidimensional slices putting them into different bracket pairs
> like `a[1..2][3..4]`. It was hell.

I worked around that in my own multidimensional array implementation by
using a helper struct for encapsulating a range of indices:

	struct IndexRange {
		int start, end;
	}

Since opSlice couldn't take more than one range, I decided to overload
opIndex instead, so you could write:

	auto arr = Array!(3,int)([3, 4, 5]);

	// Index single element:
	arr[0, 1, 2] = ...;

	// Subdimensional slicing:
	auto brr = arr[IndexRange(0,2), 2, IndexRange(3,4)];

	// opDollar works correctly in all dimensions, even inside
	// IndexRange!
	auto crr = arr[$-1, IndexRange(0,$), IndexRange($-2, $-1)];
	// equivalent to arr[2, IndexRange(0,4), IndexRange(3,4)]

This is implemented by using variadics since int differs from
IndexRange:

	auto opIndex(A...)(A args) {
		foreach (i; args) {
			static if (is(typeof(i) == int)) {
				...
			} else static if (is(typeof(i) == IndexRange)) {
				...
			}
		}
		return result;
	}

No temporaries are created, unlike the chained opSlice approach, you can
make subdimensional slices in a single opIndex call. (The disadvantage,
though, is that it has template bloat if you do a lot of this.)


> There is pull request [2] that adds everything but strides.
> My code relies on it. It still compiles even without these features
> but one have to construct slices manually: `a[Slice(1,2),
> Slice(3,4)]`.
> 
> Let's hope work on linear algebra support will legitimate the new
> features.
[...]

Yeah, I'm hoping for solid support for linear algebra eventually. D's
generics makes it possible to write generic linear algebra algorithms
that don't depend on exact matrix representation, which in C++ makes it
a pain to interface two different algebra libraries together (you have
to keep converting between representations).


T

-- 
If you look at a thing nine hundred and ninety-nine times, you are
perfectly safe; if you look at it the thousandth time, you are in
frightful danger of seeing it for the first time. -- G. K. Chesterton