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