-
Notifications
You must be signed in to change notification settings - Fork 8
/
ek_struct.m
141 lines (118 loc) · 4.44 KB
/
ek_struct.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
function [S, ST] = ek_struct(A, cholesky)
%BUILD_RK_STRUCT Build a struct for building rational Krylov subspaces.
%
if issparse(A)
if ~exist('cholesky', 'var')
cholesky = issymmetric(A);
end
if cholesky
[R, g, pA] = chol(A, 'vector');
if (g ~= 0)
% warning('Matrix is not posdef, resorting to LU');
[S, ST] = ek_struct(A, false);
return;
end
qA(pA) = 1:size(A,1);
S = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, R', R, pA, qA, x), ...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', normest(A, 1e-2));
ST = S;
else
[LA, UA, pA, qA] = lu(A, 'vector');
iqA(qA) = 1:size(A,1);
ipA(pA) = 1:size(A,1);
S = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, LA, UA, pA, iqA, x), ...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', normest(A, 1e-2));
if nargout >= 2
ST = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, UA', LA', qA, ipA, x), ...
'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
'isreal', S.isreal, ...
'nrm', S.nrm);
end
end
elseif isa(A, 'hodlr')
[LA,UA] = lu(A);
S = struct(...
... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
'solve', @(nu, mu, x) lu_solve(nu, mu, LA, UA, x),...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', norm(A));
if nargout >= 2
ST = struct(...
...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
'solve', @(nu, mu, x) lu_solve(nu, mu, UA', LA', x), ...
'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
'isreal', S.isreal, ...
'nrm', S.nrm);
end
elseif isa(A, 'hss')
ULV = ulv(A);
S = struct(...
... % 'solve', @(nu, mu, x) (nu * A - mu * eye(size(A), 'like', A)) \ x, ...
'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULV, x),...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', norm(A));
if nargout >= 2
ULVT = ulv(A');
ST = struct(...
...% 'solve', @(nu, mu, x) (nu * A' - mu * eye(size(A), 'like', A)) \ x, ...
'solve', @(nu, mu, x) ulv_solve_wrapper(nu, mu, ULVT, x), ...
'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
'isreal', S.isreal, ...
'nrm', S.nrm);
end
else
% In this case we are probably dealing with a dense matrix, so we try
% Cholesky first in the symmetric case, and otherwise resort to a
% standard LU factorization
if exist('cholesky', 'var') && cholesky
[R, g] = chol(A);
pA = (1 : length(A));
if (g ~= 0)
% warning('Matrix is not posdef, resorting to LU');
[S, ST] = ek_struct(A, false);
return;
end
S = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, R', R, pA', pA, x), ...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', normest(A, 1e-2));
ST = S;
else
% We use this syntax to be able to call sparse_solve even if this
% is the case of a generic matrix.
[LA, UA, pA] = lu(A, 'vector');
qA = 1 : length(A);
S = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, LA, UA, pA, qA, x), ...
'multiply', @(rho, eta, x) rho * (A * x) - eta * x, ...
'isreal', isreal(A), ...
'nrm', normest(A, 1e-2));
if nargout >= 2
ST = struct(...
'solve', @(nu, mu, x) sparse_solve(nu, mu, UA', LA', pA', qA', x), ...
'multiply', @(rho, eta, x) rho * (A' * x) - eta * x, ...
'isreal', S.isreal, ...
'nrm', S.nrm);
end
end
end
try
% Not all datatypes implement issymmetric at the moment
S.issymmetric = issymmetric(A);
catch
S.issymmetric = false;
end
if nargout > 1
ST.issymmetric = S.issymmetric;
end
end