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
|
function test26
%TEST26 test cs_dmsol and cs_dmspy
%
% Example:
% test26
% See also: testall
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
clear functions
randn ('state', 0) ;
rand ('state', 0) ;
clf
ntrials = 1000 ;
e1 = zeros (ntrials,1) ;
e2 = zeros (ntrials,1) ;
for trials = 1:ntrials
m = fix (100 * rand (1)) ;
n = fix (100 * rand (1)) ;
% d = 0.1 * rand (1) ;
d = rand (1) * 4 * max (m,n) / max (m*n,1) ;
A = sprandn (m,n,d) ;
% S = sprandn (m,m,d) + speye (m) ;
if (~ispc)
if (rand ( ) > .5)
A = A + 1i * sprand (A) ;
end
end
subplot (1,3,2) ; spy (A) ;
subplot (1,3,3) ; cs_dmspy (A) ;
b = rand (m,1) ;
if (~ispc)
if (rand ( ) > .5)
b = b + 1i * rand (size (b)) ;
end
end
% MATLAB cannot do A\b when A is sparse and rectangular and either
% A or b are complex
if (m ~= n & isreal (A) & ~isreal (b)) %#ok
x1 = (A\real(b)) + 1i * (A\imag(b)) ;
err1 = norm (A*x1-b) ;
elseif ((m ~= n) & ~isreal (A)) %#ok
err1 = 1 ;
else
x1 = A\b ;
err1 = norm (A*x1-b) ;
end
x2 = cs_dmsol (A,b) ;
err2 = norm (A*x2-b) ;
lerr1 = log10 (max (err1, eps)) ;
lerr2 = log10 (max (err2, eps)) ;
fprintf ('rank: %3d %3d err %6.2e %6.2e : %6.1f\n', ...
sprank(A), rank(full(A)), err1, err2, lerr1 - lerr2) ;
if (isnan (err1))
lerr1 = 10 ;
end
if (isnan (err2))
lerr2 = 10 ;
end
if (lerr2 > lerr1 + 5)
% pause
end
e1 (trials) = lerr1 ;
e2 (trials) = lerr2 ;
subplot (1,3,1) ; plot (e1, e2, 'o', [-16 10], [-16 10], 'r') ;
xlabel ('MATLAB error') ;
ylabel ('dmsol error') ;
drawnow
% pause
end
|
