Vectors and matrices

[Lars Kyllingstad]

Andrei Alexandrescu wrote:
> Lars Kyllingstad wrote:
>> I am writing a D library based some of the stuff in SLATEC, and I've 
>> come to a point where I need to decide on a way to manipulate vectors 
>> and matrices. To that end, I have some ideas and questions I would 
>> like comments on from the community.
>>
>> Ideally, I want to restrict the user as little as possible, so I'm 
>> writing heavily templated code in which one can use both 
>> library-defined vector/matrix types and built-in arrays (both static 
>> and dynamic). My reasons for this are:
>>
>>    a) Different problems may benefit from different types. Sparse 
>> matrices, dense matrices, triangular matrices, etc. can all be 
>> represented differently based on efficiency and/or memory requirements.
> 
> I use all of the above. It would be great to have them all within an 
> integrated framework.

I don't think I'm the man to provide you with that, at least not yet. I 
have next to no experience with high-performance linear algebra. This is 
a major part of the reason why I want to let people choose for 
themselves what matrix types/libraries they want to use in conjunction 
with my stuff.

Therefore I am, for now, focusing on other things. I am about halfway 
done writing D versions of the QUADPACK routines, and I have started 
work on MINPACK. It is for the latter that some basic linear algebra 
functionality is needed.

You could take a look at Bill Baxter's MultiArray library, which 
contains wrappers for several linear algebra libraries, as well as 
storage formats for different types of matrices. Also, I think, it works 
with Don's BLADE library. (Kudos to Don for that awesome name, by the 
way -- it's a lot cooler than BLAS.) Both are written in D1, though.


>>    b) I hope that, at some point, my library will be of such a quality 
>> that it may be useful to others, and in that event I will release it. 
>> Interoperability with other libraries is therefore a goal for me, and 
>> a part of this is to let the user choose other vector/matrix types 
>> than the ones provided by me.
> 
> Yes please. It would be great if you considered submitting it to Phobos.

I agree that a set of decent vector and matrix types should be put into 
phobos. But the packages I mention above are perhaps more suited for a 
dedicated scientific library, don't you think?

Right now I use D1, but I've been looking at the new phobos lately, and 
I have to say that the combination D2+phobos2 is a very enticing one. 
Writing array wrappers has never been easier than with the new "alias 
this" feature. :)

If it weren't for the lack of a 64-bit compiler for Linux I would switch 
immediately. As it is, I am seriously considering using the 32-bit one, 
even though it just feels... wrong.

Actually, I tried compiling my lib with the latest 32-bit DMD2. 
Immediately after pressing enter, memory usage went through the roof and 
my computer became completely unresponsive. It took me forever to kill 
the dmd process. I guess it has something to do with the heavy use of 
templates. Has anyone else experienced this?


>>    c) Often, for reasons of both efficiency and simplicity, it is 
>> desirable to use arrays directly.
> 
> Yah.
> 
>> My first question goes to those among you who do a lot of linear 
>> algebra in D: Do you think supporting both library  types and arrays 
>> is worth the trouble? Or should I just go with one and be done with it?
> 
> If you go templated you don't need to explicitly support built-in arrays 
> - they'll just work.

Yeah, vectors work fine, especially now that D has built-in vector 
operations. The problem is with matrices, as described in the following.


>> A user-defined matrix type would have opIndex(i,j) defined, and to 
>> retrieve elements one would write m[i,j].
> 
> Yah.
> 
>> However, the syntax for two-dimensional arrays is m[i][j], and this 
>> means I have to put a lot of static ifs around my code, in order to 
>> check the type every time I access a matrix.
> 
> No, that's not a two-dimensional array; it's an array of arrays. If you 
> want to make your lib work with arrays of arrays, you could easily build 
> a little wrapper arround it (e.g. JaggedMatrix).

OK, technically it may not be what is called a two-dimensional array. 
But it's the closest we've got, no? Or perhaps that would be a 
rectangular (static) array, but those still have the [][] indexing syntax.


>> This leads me to my second question, which is a suggestion for a 
>> language change, so I expect a lot of resistance. :)
>>
>> Would it be problematic to define m[i,j,...] to be equivalent to 
>> m[i][j][...] for built-in arrays, so that arrays and user-defined 
>> types could be used interchangeably?
>>
>> (And, importantly, are there anyone but me who think they would 
>> benefit from this?)
> 
> This wouldn't harm, but it would be a special case.

I know, and possibly an esoteric one, at that. But it would make it 
easier to create drop-in array replacements.

-Lars