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
|
function demo2 (C, sym, name)
%DEMO2: solve a linear system using Cholesky, LU, and QR, with various orderings
%
% Example:
% demo2 (C, 1, 'name of system')
% See also: cs_demo
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
clf
subplot (2,2,1) ; cspy (C) ;
title (name, 'FontSize', 16, 'Interpreter', 'none') ;
[m n] = size (C) ;
[p,q,r,s,cc,rr] = cs_dmperm (C) ;
subplot (2,2,3) ; cs_dmspy (C) ;
subplot (2,2,2) ; ccspy (C) ;
drawnow
sprnk = rr (4) - 1 ;
nb = length (r) - 1 ;
ns = sum ((r (2:nb+1) == r (1:nb)+1) & (s (2:nb+1) == s (1:nb)+1)) ;
fprintf ('blocks: %d singletons %d structural rank %d\n', nb, ns, sprnk) ;
if (sprnk ~= sprank (C))
error ('sprank mismatch!') ;
end
if (sprnk < min (m,n))
return ; % return if structurally singular
end
% the following code is not in the C version of this demo:
if (m == n)
if (sym)
try
[L,p] = cs_chol (C) ; %#ok
subplot (2,2,4) ;
cspy (L+triu(L',1)) ; title ('L+L''') ;
catch
% tol = 0.001 ;
[L,U,p,q] = cs_lu (C,0.001) ; %#ok
subplot (2,2,4) ;
cspy (L+U-speye(n)) ; title ('L+U') ;
end
else
[L,U,p,q] = cs_lu (C) ; %#ok
subplot (2,2,4) ;
cspy (L+U-speye(n)) ; title ('L+U') ;
end
else
if (m < n)
[V,beta,p,R,q] = cs_qr (C') ; %#ok
else
[V,beta,p,R,q] = cs_qr (C) ; %#ok
end
subplot (2,2,4) ;
cspy (V+R) ; title ('V+R') ;
end
drawnow
% continue with the MATLAB equivalent of the C cs_demo2 program
for order = [0 3]
if (order == 0 & m > 1000) %#ok
continue ;
end
fprintf ('QR ') ;
print_order (order) ;
b = rhs (m) ; % compute right-hand-side
tic ;
x = cs_qrsol (C, b, order) ;
fprintf ('time %8.2f ', toc) ;
print_resid (C, x, b) ;
end
if (m ~= n)
return ;
end
for order = 0:3
if (order == 0 & m > 1000) %#ok
continue ;
end
fprintf ('LU ') ;
print_order (order) ;
b = rhs (m) ; % compute right-hand-side
tic ;
x = cs_lusol (C, b, order) ;
fprintf ('time %8.2f ', toc) ;
print_resid (C, x, b) ;
end
if (sym == 0)
return ;
end
for order = 0:1
if (order == 0 & m > 1000) %#ok
continue ;
end
fprintf ('Chol ') ;
print_order (order) ;
b = rhs (m) ; % compute right-hand-side
tic ;
x = cs_cholsol (C, b, order) ;
fprintf ('time %8.2f ', toc) ;
print_resid (C, x, b) ;
end
|
