the competition's expm vs ours

[Marco Caliari]

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
-------------- next part --------------
function [E,AA] = expmdemo1(A)
% trace reduction
A(find(A==-Inf))=-realmax;
trshift = sum(diag(A))/length(A);
if (trshift > 0)
  A = A-trshift*speye(size(A));
end
% balancing
[DP,AA] = balance(A,"P");
[DS,AAA] = balance(AA,"S");
AA = AAA;
[f,e] = log2(norm(AA,'inf'));
s = max(0,e);
s = min(s,1023);
AA = AA/2^s;
% Pade approximation for exp(A)
c = [5.0000000000000000e-1,...
     1.1666666666666667e-1,...
     1.6666666666666667e-2,...
     1.6025641025641026e-3,...
     1.0683760683760684e-4,...
     4.8562548562548563e-6,...
     1.3875013875013875e-7,...
     1.9270852604185938e-9];
q = 8;
N = sparse(zeros(size(AA)));
D = N;
p = 1;
for k = q:-1:1
  N = N+c(k)*speye(size(AA));
  N = AA*N;
  D = D+p*c(k)*speye(size(AA));
  D = AA*D;
  p = -p;
end
N = N+speye(size(AA));
D = D+speye(size(AA));
E = full(D\N);
% Undo scaling by repeated squaring
for k = 1:s 
  E = E*E;
end
% inverse balancing
E = DS*E/DS;
E = DP*E*DP';
% inverse trace reduction
if (trshift >0)
  E = E*exp(trshift);
end