ch-ch-update: series, closed-form series, and strides

[Andrei Alexandrescu]

Bill Baxter wrote:
> On Sat, Jan 31, 2009 at 12:12 PM, Andrei Alexandrescu
> <SeeWebsiteForEmail at erdani.org> wrote:
>> I've updated my code and documentation to include series (as in math) in the
>> form of infinite ranges. Also series in closed form (given n can compute the
>> nth value without iterating) are supported as random-access ranges.
>>
>> Also Stride is provided. The Matrix container (speaking of scientific
>> computing with D!) will support various representational choices, most
>> importantly the ones endorsed by high-performance libraries. For Matrix,
>> Stride is an important component as I'm sure anyone who's ever written a
>> matrix knows.
> 
> This sounds interesting!  So you'll have diagonal, banded, upper/lower
> triangular, symmetric, strided and dense matrices?  How about 2D
> slices?

Sparse matrices too. For dense block matrices I'm thinking of a type 
like this:

alias BlockMatrix!(float, 3, [0, 2, 1]) M3;

This defines a block rectangular matrix with three dimensions laid out 
in the order specified: In this case, M3 is laid out as a collection of 
column-order submatrices of 2 dimensions. (I deliberately chose a 
slightly odd example.) The last parameter of BlockMatrix can be any 
permutation of 0, 1, and 2, thus dictating how the dimensions are laid 
out. For an obvious example, here's a Fortran-compatible block matrix 
with column-major layout:

alias BlockMatrix!(double, 2, [1, 0]) FortranMatrix;

So it's pretty easy to choose any layout for any matrix with any number 
of dimensions. This matrix defines ranges that crawl it, called 
hyperplanes. There is a "natural" hyperplane that doesn't need any 
stride, and that's the cut along the dimension that varies slowest. That 
cut is the most efficient because it really has the layout of a 
(littler) BlockMatrix itself. Continuing on the M3 example, a "row" in 
that matrix has type:

alias BlockMatrix!(float, 2, [1, 0]) M3Row;

which is obtained by chopping the first element in the array [0, 2, 1] 
into [2, 1], and then decrementing what's left into [1, 2]. So if you have:

M3 m;
auto subM = m[3];

you get that smaller BlockMatrix with all amenities its larger cousin 
offered. But if you want a cut along a different dimension, a strided 
range will be returned.

Jagged, banded, diagonal, fixed-size, and sparse layouts come naturally 
and will be hopefully supported in a mixable-and-matchable form (e.g. I 
want a block matrix of 3x3 matrices). I want to focus on storage for 
now, put it in the standard library, to then allow great scientific 
programmers like you guys to use those storages as a common data format 
on top of which cool algorithms can be implemented.


Andrei