-
Notifications
You must be signed in to change notification settings - Fork 7
/
logm2.m
54 lines (49 loc) · 1.33 KB
/
logm2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
function [Y, warning] = logm2(A, pert)
%LOGM2 LOGM2(X) is the matrix natural logarithm of X . Complex
% results are produced if X has nonpositive eigenvalues.
% LOGM2 can be thought of as being computed using eigenvalues
% D and eigenvectors V, such that if [V,D] = EIG(X) then
% LOGM2(X) = V*log(D)/V.
% Copyright (c) 1986-93 by the MathWorks, Inc.
warning = 0;
if ~length(A), Y = []; return; end
[Q,T] = schur(A);
[Q,T] = rsf2csf(Q,T);
dT = diag(T);
F = log(diag(T));
% Set log of zero to large negative number
find_vec = find(~finite(F));
F(find_vec) = -ones(length(find_vec),1).*(1/eps);
F = diag(F);
tol = eps*norm(T,1);
n = max(size(A));
for p = 1:n-1
for i = 1:n-p
j = i+p;
s = T(i,j)*(F(j,j)-F(i,i));
if p > 1
k = i+1:j-1;
s = s + T(i,k)*F(k,j) - F(i,k)*T(k,j);
end
d = T(j,j) - T(i,i);
if abs(d) <= tol,
warning = 1;
d = tol;
end
F(i,j) = s/d;
end
end
Y = Q*F*Q';
% If diagonal elements are the same then check accuracy against expm.
% If accuracy is poor then perturb matrix and call logm2 again.
if (warning)
tol = tol + eps;
if norm(expm(Y) - A, 'inf') > 1e5*tol
if nargin ==1,
pert = tol;
else
pert = pert * 100
end
Y = logm2(A + pert*randn(n,n), pert);
end
end