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