I've put the beginnings of my matrix-vector library up on github (between swizzling, generalising to N dimensions and making use of static foreach this project quickly turned from a cleanup to a rewrite). You can find it here.

http://github.com/gcharnock/phoboslinalgebra

Some highlights so far:

//Dimension and field are templatable parameters. There is no hard limit on N, but obviously everything is set up for N being small.
alias Vector!(float,3) vector_t;
alias Matrix!(float,3) matrix_t;

//Very natural initialisation
auto v=vector_t(1,2,3);
auto m=matrix_t(1,2,0,0,1,0,0,3,1);

//Single element access does not require opDispatch
writeln(v.x);

//Swizzling supported via opDispatch
writeln(v.yzx); //2 3 1
writeln(v.xxyyzx); //1 1 2 2 3 1

//No matter how big N, elements always have a canonical name (probably more important for matrices where the letters of the alphabet run out faster although I've only done this for vectors so far)
writeln(v.x2); //3
writeln(v.x0x1x2); //1 2 3, zero based indexing because D inherits C

//Some Math operations. Static foreach is used to unroll all loops. Unfortunately I ran into some problems putting the operators outside the struct. I'm not sure if this was me or the compiler but I'll need another hacking session to work out what's going on.

float s;
v*s    v/s    v*=s
v/=s   v+v    v-v
v*v    m+m    m-m
m*m    m*v

A question I'd like input on would what should the semantics of the Transpose(), Transposed() function. I can think of some alternatives:

Transpose:
1) Naive solution, copies N^^2-N scalars.
2) Set a bool. Switch between two element iteration schemes for nearly every operation (with suitable use of templates this might not be as painful as it sounds.)

Transposed:
1) makes a new matrix
2) creates a transposed image which is forever linked to the original matrix: A=B.Transpose() => Changes to A affect B
3) creates a transposed image which is copy-on-write linked to the original matrix: A=B.Transpose() => Changes to A do not affect B

On Sun, 25 Apr 2010 23:17:10 -0300, Gareth Charnock <gareth.tpc@gmail.com> wrote:
> I've put the beginnings of my matrix-vector library up on github (between swizzling, generalising to N dimensions and making use of static foreach this project quickly turned from a cleanup to a rewrite). You can find it here.
>
> http://github.com/gcharnock/phoboslinalgebra
>
> Some highlights so far:
>
> //Dimension and field are templatable parameters. There is no hard limit on N, but obviously everything is set up for N being small.
> alias Vector!(float,3) vector_t;
> alias Matrix!(float,3) matrix_t;
>
> //Very natural initialisation
> auto v=vector_t(1,2,3);
> auto m=matrix_t(1,2,0,0,1,0,0,3,1);
>
> //Single element access does not require opDispatch
> writeln(v.x);
>
> //Swizzling supported via opDispatch
> writeln(v.yzx); //2 3 1
> writeln(v.xxyyzx); //1 1 2 2 3 1
>
> //No matter how big N, elements always have a canonical name (probably more important for matrices where the letters of the alphabet run out faster although I've only done this for vectors so far)
> writeln(v.x2); //3
> writeln(v.x0x1x2); //1 2 3, zero based indexing because D inherits C
>
> //Some Math operations. Static foreach is used to unroll all loops. Unfortunately I ran into some problems putting the operators outside the struct. I'm not sure if this was me or the compiler but I'll need another hacking session to work out what's going on.
>
> float s;
> v*s    v/s    v*=s
> v/=s   v+v    v-v
> v*v    m+m    m-m
> m*m    m*v
>
> A question I'd like input on would what should the semantics of the Transpose(), Transposed() function. I can think of some alternatives:
>
> Transpose:
> 1) Naive solution, copies N^^2-N scalars.
> 2) Set a bool. Switch between two element iteration schemes for nearly every operation (with suitable use of templates this might not be as painful as it sounds.)
>
> Transposed:
> 1) makes a new matrix
> 2) creates a transposed image which is forever linked to the original matrix: A=B.Transpose() => Changes to A affect B
> 3) creates a transposed image which is copy-on-write linked to the original matrix: A=B.Transpose() => Changes to A do not affect B

Some quick comments:
- The version I'm seeing on github doesn't seem to have all the features you're referencing (i.e. v*v). Why are some ops limited to */ and other +-?
- Static foreach isn't a good way to do loop unrolling (at least historically).
- In terms of naming, I prefer float3 instead of vector3f.
- Return type and argument types need to be dynamically deduced: i.e. int3 * float3 should work.
- length should mean the number of vector elements. I'd recommend norm or L2 instead.
- For the matrix, Field[D][D] raw; will let you get out of manual indexing issues.
- I've never seen matrix swizzling and don't really know what it would be used for.
- You should be able to set matrix layout to either C or Fortran style
- Re: Transpose: this should be an in place swap. Given the value-type semantics of small matrices, make the copies.
- Re: Transposed: this function should be named T (or have an equivalent alias) and return a reference to the current matrix casted to its C or Fortran layout counterpart.

In general, this feels like a decent start, although it's still very alpha and feature incomplete. Keep up the good work.

Thanks. To quickly answer this:

> - The version I'm seeing on github doesn't seem to have all the features
> you're referencing (i.e. v*v). Why are some ops limited to */ and other +-?

It was quite late (3am) when I typed up that email. I'm sorry if I got the ops wrong. v*v was actually rejected in favour of dot(v,v). Unless I've made another mistake the implemented operators are:

v*s    v/s    v*=s
v/=s   v+v    v-v
v+=v   v-=v   m+m
m-m    m*m    m*v

(where v is a vector, m a matrix, s a scalar)

The reason the feature set is very limited is I wanted to get feedback before implementing a large set of features and then having to change things. The idea was to implement just enough so that it would be clear where I was going. For example, I wasn't aware that compile time foreach was a bad way to unroll loops (perhaps I should build a mixin string?).

Gareth Charnock:

>http://github.com/gcharnock/phoboslinalgebra

I think that a single module is better. If seen fitting, some functionality can be moved inside other already present modules of Phobos.
The module will need a good amount of unit tests.

In the code I see no point in calling functions like:
CTFENextToken
Just call them:
nextToken
Note: functions, template functions, and ctfe functions, start with a lower case letter (and templates generally with with an upper case).

Don't call this module and its contents matrix or vector or Vector, etc. Call them SmallVector, or SmallVec, SmallMat, ecc, because this lib is not designed to be a general purpose matrix or vector lib, it's designed for small number of items.
You can call this module std.tinyarray :-)


> alias Matrix!(float,3) matrix_t;

To define the sizes of a matrix you need two numbers:
alias Matrix!(float, 3, 5) Matrix35;


> //Very natural initialisation
> auto m=matrix_t(1,2,0,0,1,0,0,3,1);

A 1D initialization of a matrix is not so natural :-)


> //No matter how big N, elements always have a canonical name (probably more important for matrices where the letters of the alphabet run out faster although I've only done this for vectors so far)

What's the point of giving names to matrix elements?


> 2) Set a bool. Switch between two element iteration schemes for nearly every operation (with suitable use of templates this might not be as painful as it sounds.)

The point of a transposition is sometimes to actually have data in a different position. With cache effects the position of data can be quite important (a matrix multiplication can be 2 or 3 times faster if you transpose. This is true for larger matrices).


> Transposed:
> 1) makes a new matrix
> 2) creates a transposed image which is forever linked to the original
> matrix: A=B.Transpose() => Changes to A affect B

The semantics of transpose/transposed can be like in Python: the transpose works in-place, transposed creates a new matrix.

Methods and properties like transpose have to start with a lowercase: auto A = B.transposed;

Bye,
bearophile

You can do performance benchmarks compared to the tinyvector module here, for D1:
http://www.fantascienza.net/leonardo/so/libs_d.zip
It comes from a performance tuning for DMD.

Bye,
bearophile

> The semantics of transpose/transposed can be like in Python: the transpose works in-place, transposed creates a new matrix.

That's what I do, with normalize/normalized too.

I'm not sure how your static foreach is actually static.

I think static foreach is a very good idea to minimize speed difference between debug and release versions (i know one shouldn't ship debug version, but sometimes you have to)...

> //Single element access does not require opDispatch
> writeln(v.x);
> 
> //Swizzling supported via opDispatch
> writeln(v.yzx); //2 3 1
> writeln(v.xxyyzx); //1 1 2 2 3 1


Can we get assignment through swizzle ? Like:
v.xy = vector2f(1.f, 2.f);

(but one should disallow: v.xx = vector2f(1.f, 2.f); )

Also, I like the naming, even though I would probably alias it anyway.

Existing solutions: vector3f, vec3f, smallVector3f, float3...

#ponce:

> I'm not sure how your static foreach is actually static.

It's a foreach on:

template TypeNuple(T, size_t n) {
    static if(n == 0) {
        alias TypeTuple!() TypeNuple;
    }
    else {
        alias TypeTuple!(T,TypeNuple!(T, n-1)) TypeNuple;
    }
}

So with the current D it's a static foreach.

Bye,
bearophile

I need a better name for ArgList. It originally made sense because because I assumed I be using the tuple like this:
this(ArgList argList)

But then I ended up using it in a load of static foreaches.

I think it would be useful to actually be able to say "static foreach" so that readers (and writers) of the code know for sure what's going on. I seem to remember you starting a thread on that very matter.

bearophile wrote:
> #ponce:
> 
>> I'm not sure how your static foreach is actually static.
> 
> It's a foreach on:
> 
> template TypeNuple(T, size_t n) {
>     static if(n == 0) {
>         alias TypeTuple!() TypeNuple;
>     }
>     else {
>         alias TypeTuple!(T,TypeNuple!(T, n-1)) TypeNuple;
>     }
> }
> 
> So with the current D it's a static foreach.
> 
> Bye,
> bearophile

Gareth Charnock:

> I think it would be useful to actually be able to say "static foreach" so that readers (and writers) of the code know for sure what's going on. I seem to remember you starting a thread on that very matter.

http://d.puremagic.com/issues/show_bug.cgi?id=4085

Bye,
bearophile