-
Notifications
You must be signed in to change notification settings - Fork 4
/
cqtadi.m
46 lines (31 loc) · 892 Bytes
/
cqtadi.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
function X = cqtadi(A, B, C, alpha, beta, varargin)
%CQTADI CQT-version of the ADI iteration
p = inputParser;
addParameter(p, 'tol', cqtoption('threshold'));
addParameter(p, 'maxit', 1000);
addParameter(p, 'debug', false);
parse(p, varargin{:});
tol = max(p.Results.tol, 10 * cqtoption('threshold'));
maxit = p.Results.maxit;
debug = p.Results.debug;
it = 1;
c = 1;
X = cqt([], []);
I = cqt(1, 1);
res = inf;
while it < maxit
Xold = X;
X = - (A - beta(c) * I) \ ( X * (B + beta(c) * I) + C );
X = - ((A - alpha(c) * I) * X + C) / (B + alpha(c) * I);
oldres = res;
res = norm(X - Xold, inf); % / norm(X, inf);
if debug
fprintf('ADI :: Iteration %d, res = %e, %e\n', it, res, norm(A*X + X*B + C, inf));
end
if res < tol % || oldres < res
break;
end
c = 1 + mod(c, length(alpha));
it = it + 1;
end
end