the competition's expm vs ours

[Marco Caliari]

Dear Jaroslav,

I can't fully test because of the balance change, but I think it is OK. I 
would only remove line 112

  a(a == -Inf) = -realmax;

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.

Best regards,

Marco

On Tue, 9 Dec 2008, Jaroslav Hajek wrote:

> On Sat, Nov 22, 2008 at 9:51 AM, Marco Caliari <marco.caliari at univr.it> wrote:
>> Dear all,
>>
>> sorry for the delay. I sent an email but probabily it went lost. In Octave
>> expm there are first some "reductions" (trace reduction and balancing) and
>> then Pade' approximation (8,8) with scaling and squaring.
>> On the other hand, Matlab does not make any reductions and uses Pade'
>> approximation (6,6) with scaling and squaring.
>>
>> I have an m function implemeting exactly what Octave does with expm (like
>> expmdemo1.m for Matlab). I can send it out of the list.
>>
>> Best regards,
>>
>> Marco
>
> Hi all,
>
> I've just written an m-file version of expm based on Marco Caliari's
> demo above, and compared
> with the built-in expm using recent tip (note that balance was
> modified to be Matlab compatible).
> Attached.
>
> A quick test on a 1000x1000 random matrix shows 5.5 seconds for the
> built-in version, 3.8 seconds for the m-file version.
> Anyone has objections against replacing?
>
> I guess an m-file version is going to be much easier to play with.
>
> regards
>
> -- 
> RNDr. Jaroslav Hajek
> computing expert
> Aeronautical Research and Test Institute (VZLU)
> Prague, Czech Republic
> url: www.highegg.matfyz.cz
>