Re: the competition's expm vs ours

[Jaroslav Hajek]

On Tue, Dec 9, 2008 at 5:54 PM, Marco Caliari <address@hidden> wrote:
> Dear Jaroslav,
>
> I can't fully test because of the balance change, but I think it is OK.

It's like a single call to xGEBAL. I was not sure whether there's a
special reason to call balance twice, which can, in rare cases, yield
slightly different results.

> I
> would only remove line 112
>
>  a(a == -Inf) = -realmax;
>

OK.

> It was just a tricky idea to manage the -Inf case (see
> http://www.nabble.com/Octave-hangup-for-expm-with-non-finite-arguments-tt14914598.html#a14960331).
> The trick with a2 is exactly what I tried some days ago to improve the code.

Yeah, that means doing it with 7 matrix multiplications instead of the
former 16, which is a win. In fact, with this optimization more time
seems to be spent in the repeated squaring. I wonder whether there's
room to optimize further, but I don't see anything left. Maybe matlab
has an even better algorithm?

cheers

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz