Re: [boost] different matrix library?

Subject: Re: [boost] different matrix library?
From: Edward Grace (ej.grace_at_[hidden])
Date: 2009-08-17 10:18:15

>> I don't quite understand why you wouldn't want a DSL like MATLAB?
>> After all, MATLAB remember is essentially FORTRAN with different
>> horse shoes, that's why so many scientific / engineering types pick
>> it up so fast - there is very little that's new to learn.
>
> Not saying that I wouldn't want it, saying that perhaps there are  
> some gems
> in other languages, too.

Ok, I misunderstood - that I agree with.  I like 'breeding' ;-)  If  
you take the best of breeds the animal you can end up with can be so  
much fitter than the originals.  Of course, along the way you might  
spawn a few freaks too!

As I've mentioned to Joel, I don't think a blind carbon copy of  
MATLAB/FORTRAN-like operation is ultimately what should be done.   
Learning from the weaknesses and addressing those (e.g. my banging on  
about tensors) has to be the way forward.

> I do not have a Fortran background, but I have to
> admit that I have used R rather extensively. R/S stuff like
>
> A[row(A)==col(A)] <- 2

Hmm, intriguing.

> is easy to use, but does not really match with something like
>
> std::fill( diag(A).begin(), diag(A).end(), 2 );

How about simply

   A = 2*identity(2);

>>   alpha(range(2,10))=alpha(range(1,9));
>>
>> or some-such?
>
> Not sure, but given your example, something like
>
> std::rotate( alpha.begin(), alpha.begin()+1, alpha.end() );

Too obscure...  Maybe that's what it could map too, but I don't think  
it's clear.  Similarly if you had some other requirement you would  
end up calling something other than 'rotate' - having a whole bunch  
of different functions to call based on a subtle indexing change  
can't be good!

> :-). I agree the Matlab code is more readable in this case.

It's actually FORTRAN - (I'm teasing of course - if there were ; s at  
the end it would have been MATLAB) ;-)

>> How about,
>>
>>   vec=mat(:,4)
>>
>> (Again, is that MATLAB or is it FORTRAN?)
>
> Are you copying or keeping a reference?

It could be either.  I'm willing to bet that in FORTRAN it will -  
behind the scenes - pick the best one contingent on what happens next  
(FORTRAN is always pass by reference for example). Why should the  
user care as long as they can do

  vec(1)

and get the correct value?

> vec = row(mat,4)

Very uBlas (nothing wrong with that). It gets tricker with > 3  
dimensions though!

>> Can you give some examples?
>>
>
> E.g., I would prefer
>
> x = solve( A, b )
>
> over x = A\b;

Sure.

>> get used just by hobbyists who enjoy the process rather than the
>> result. (I'm guilty of that myself - hence I'm here).
>
> I'm not after that. I'm just saying that there might be people who  
> actually
> do stuff in other packages than Matlab (like R et al.).

Ok. I agree - 'breeding'!

> What I also wanted to add is that many (real-life) examples of  
> algorithms do
> not contain that much linear algebra.

Usually a single call to 'solve' (in your parlance), in my experience  
the rest of the code is setting up the problem!

> I've seen a great deal of .m files
> which have about 20% linear algebra. About 80% is conditional  
> program flow
> and bookkeeping stuff. This is where C++ has a clear edge. If we're  
> going to
> mix them, it should look natural.

Sure.  Most of my above discussion is centred around setting up the  
problem.  For example you need to build some sort of complicated  
operator (Matrix) to then apply to a state vector.

The tricky part is elegantly expressing how to build that - being  
able to easily refer to subdomains, etc is essential.

A concrete example of what I mean is a little function I wrote  
'fftnpad' here,

http://www.boostpro.com/vault/index.php? 
direction=0&order=&directory=linalg&

this takes an N-dimensional data structure and pads it assuming FFT  
ordering -- inserting blank space in the middle.  In 1D this is  
trivial, 2D -- ok that's simple enough, how about 4D? (Yes I do  
occasionally do that).

It's conceptually not very difficult but describing the regions you  
want can be tricky.  The secret is the

   B(tgt{:}) = A(src{:})

at the end which makes life very easy. How would one code this in C++  
for example?  Give it a bash - it's trickier than it looks (I'm happy  
for tips to do this better by the way).

Perhaps, as a general suggestion, it's worthwhile building a  
repository of 'gymnastics' such as this in various different  
languages.  That may prove invaluable to people in the future  
'breeding' process.

-ed
------------------------------------------------
"No more boom and bust." -- Dr. J. G. Brown, 1997