-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfunmv.m
190 lines (176 loc) · 6.22 KB
/
funmv.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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
function [C,S,s,m,mv,mvd,unA] = funmv(t,A,b,flag,M,prec,shift,bal,full_term)
%FUNMV trigonometric and hyperbolic matrix functions times vectors
% [C,S,s,m,mv,mvd] = FUNMV(t,A,B,FLAG,[],PREC) computes f(t*A)*B or
% f(t*SQRTM(A))B, where f is a trigonometric or hyperbolic matrix function
% and SQRTM denotes any matrix square root, without
% explicitly forming f(t*A) or SQRTM(A). PREC is the required accuracy,
% 'double', 'single' or 'half', and defaults to CLASS(A).
% A total of mv products with A or A^* are used, of which mvd are
% for norm estimation.
% The full syntax is
% [C,S,s,m,mv,mvd,unA] = funmv(t,A,b,flag,M,prec,shift,bal,full_term).
% unA = 1 if the alpha_p were used instead of norm(A).
% If repeated invocation of FUNMV is required for several values of t
% or B, it is recommended to provide M as an external parameter as
% M = SELECT_TAYLOR_DEGREE_TRIG(A,b,flag,m_max,p_max,prec,shift,bal,true).
% This also allows choosing different m_max and p_max.
% The outputs C and S depend on the input FLAG as follows:
% FLAG = 'cos.sin' , C = cos(tA)B and S = sin(tA)B
% FLAG = 'cosh.sinh' , C = cosh(tA)B and S = sinh(tA)B
% FLAG = 'cos.sinc' , C = cos(tA)B and S = sinc(tA)B
% FLAG = 'cosh.sinch' , C = cosh(tA)B and S = sinch(tA)B
% FLAG = 'cos.sinc.sqrt' , C = cos(t*sqrt(A))B and S = sinc(t*sqrt(A))B
% FLAG = 'cosh.sinch.sqrt', C = cosh(t*sqrt(A))B and S = sinch(t*sqrt(A))B
% The parameter SHIFT is optinal only if FLAG = 'cos.sin'.
% Example:
% To evaluate a combination of the form
% cos(t*sqrt(A))*y0 + t*sinc(t*sqrt(A))*y1,
% it is C(:,1) + t*S(:,2), where [C,S] = funmv(t,A,[y0,y1],'cos.sinc.sqrt')
%
% Reference: A. H. Al-Mohy, A truncated Taylor series algorithm for
% computing the action of trigonometric and hyperbolic matrix
% functions, SIAM J. Sci. Comp., 40(3), A1696-A1713.
% Awad H. Al-Mohy, January 18, 2018.
if nargin < 4
error('funmv:NoInput', 'Not enough input parameters to funmv.')
end
if ~ismember(flag, {'cos.sin', 'cosh.sinh', 'cos.sinc'...
'cosh.sinch', 'cos.sinc.sqrt', 'cosh.sinch.sqrt'})
error('funmv:NoInput', 'Choose a correct input of the parameter FLAG.')
end
if nargin < 7 || isempty(shift), shift = true; end
flag1 = isequal( flag, 'cos.sin');
flag2 = isequal( flag, 'cosh.sinh');
if flag1 , flag = 1; sign = 1; end
if flag2 , flag = 1; sign = 0; end
if isequal( flag, 'cos.sinc') , flag = 1; sign = 1; end
if isequal( flag, 'cosh.sinch') , flag = 1; sign = 0; end
if isequal( flag, 'cos.sinc.sqrt') , flag = 0; sign = 1; end
if isequal( flag, 'cosh.sinch.sqrt'), flag = 0; sign = 0; end
if ~(flag1 || flag2), shift = false; end
if nargin < 9 || isempty(full_term), full_term = false; end
if nargin < 8 || isempty(bal), bal = false; end
if bal
[D,B] = balance(A);
if norm(B,1) < norm(A,1), A = B; b = D\b; else bal = false; end
end
n = length(A); n0 = size(b,2);
pp = n0;
if flag, pp = 2*n0; end % the number of matrix-vector prod., A(AB).
if shift && (flag1 || flag2)
mu = trace(A)/n;
mu = full(mu); % The code could be slower without the full!
tmu = t*mu;
A = A - mu*speye(n);
end
if nargin < 6 || isempty(prec), prec = class(A); end
t2 = t;
if ~flag, t2 = t^2; end
if nargin < 5 || isempty(M)
tt = 1;
[M,mvd,~,unA] = select_taylor_deg_trig(t2*A,b,flag,[],[],prec,false,false,false);
mv = mvd;
else
tt = t; mv = 0; mvd = 0; unA = 0;
end
switch prec
case 'double', tol = 2^(-53);
case 'single', tol = 2^(-24);
case 'half', tol = 2^(-10);
end
s = 1;
if t == 0
m = 0;
else
[m_max,p] = size(M);
S = diag(1:m_max);
C = ( (ceil(abs(tt)*M))'*S );
C (C == 0) = inf;
if p > 1
[cost, m] = min(min(C)); % cost is the overall cost.
else
[cost, m] = min(C); % when C is one column. Happens if p_max = 2.
end
if cost == inf; cost = 0; end
s = max(cost/m,1);
end
undo_inside = false; undo_outside = false;
% undo shifting inside or outside the loop
if shift
if flag1 && ~isreal(tmu)
cosmu = cos(tmu/s); sinmu = sin(tmu/s); undo_inside = true;
elseif flag1 && isreal(tmu) && abs(tmu)> 0
cosmu = cos(tmu); sinmu = sin(tmu); undo_outside = true;
elseif flag2 && abs(real(tmu)) > 0
cosmu = cosh(tmu/s); sinmu = sinh(tmu/s); undo_inside = true;
elseif flag2 && ~real(tmu) && abs(tmu)> 0
cosmu = cosh(tmu); sinmu = sinh(tmu); undo_outside = true;
end
end
mods = mod(s,2);
C0 = zeros(n,n0);
if mods, C0 = b/2; end
S = C0;
C1 = b;
for i = 1:s+1
if i == s+1
S = 2*S; C1 = S;
end
V = C1;
if undo_inside, Z = C1; end
b = C1;
c1 = norm(b,inf);
for k = 1:m
even = 2*k;
if i <= s
odd = even-1; q = 1/(even+1);
else % when i = s+1, compute Taylor poly. of sinc
odd = even+1; q = odd;
end
if flag, b = A*b; end
b = (A*b)*((t/s)^2/(even*odd));
mv = mv + pp;
V = V + ((-1)^(k*sign))*b;
if undo_inside, Z = Z + (((-1)^(k*sign))*q)*b; end
c2 = norm(b,inf);
if ~full_term
if c1 + c2 <= tol*norm(V,inf)
break
end
c1 = c2;
end
end
if undo_inside
if i <= s
V = V*cosmu + A*(Z*(((-1)^sign)*t*sinmu/s));
mv = mv + n0;
else
V = A*(V*(t*cosmu/s)) + Z*sinmu;
mv = mv + n0;
end
end
if i == 1
C2 = V;
elseif i <= s
C2 = 2*V - C0;
end
if i <= s-1 && xor( mods,mod(i,2) )
S = S + C2; % sum of C_i for even i's if s is odd and vice versa
end
C0 = C1; C1 = C2;
end
C = C2;
if undo_inside
S = V;
elseif flag1 || flag2
S = A*(V*(t/s));
mv = mv + n0;
else
S = V/s; % sinc(tA)B (or sinch(tA)B ) if flag = 1 and
% sinc(t*sqrt(A))B (or sinch(t*sqrt(A))B) if flag = 0
end
if undo_outside
C = cosmu*C + (((-1)^sign)*sinmu)*S;
S = sinmu*C2 + cosmu*S;
end
if bal, C = D*C; S = D*S; end