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
|
function test6
%TEST6 test cs_reach, cs_reachr, cs_lsolve, cs_usolve
%
% Example:
% test6
% See also: testall
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
rand ('state', 0)
maxerr = 0 ;
clf
for trial = 1:201
n = fix (100 * rand (1)) ;
d = 0.1 * rand (1) ;
L = tril (sprandn (n,n,d),-1) + sprand (speye (n)) ;
b = sprandn (n,1,d) ;
if (~ispc)
if (rand ( ) > .5)
L = L + 1i*sprand(L) ;
end
if (rand ( ) > .5)
b = b + 1i*sprand(b) ;
end
end
for uplo = 0:1
if (uplo == 1)
% solve Ux=b instead ;
L = L' ;
end
x = L\b ;
sr = 1 + cs_reachr (L,b) ;
sz = 1 + cs_reachr (L,b) ;
check_if_same (sr,sz) ;
s2 = 1 + cs_reach (L,b) ;
try
if (uplo == 0)
x3 = cs_lsolve (L,b) ;
else
x3 = cs_usolve (L,b) ;
end
catch
if (isreal (L) & isreal (b)) %#ok
lasterr
error ('!') ;
end
% punt: sparse(L)\sparse(b) not handled by cs_lsolve or cs_usolve
x3 = L\b ;
end
% x3 is NOT returned in sorted order, so it is not a valid MATLAB
% sparse matrix. It might also have explicit zeros. Double-transpose
% it to sort it, and multiply one to remove explicit zeros.
x3 = 1 * (x3')' ;
spy ([L b x x3])
drawnow
s = sort (sr) ;
[i j xx] = find (x) ; %#ok
[i3 j3 xx3] = find (x3) ; %#ok
if (isempty (i))
if (~isempty (s))
i %#ok
s %#ok
error ('!') ;
end
elseif (any (s ~= i))
i %#ok
s %#ok
error ('!') ;
end
if (isempty (i3))
if (~isempty (s))
i3 %#ok
s %#ok
error ('!') ;
end
elseif (any (s ~= sort (i3)))
s %#ok
i3 %#ok
error ('!') ;
end
if (any (s2 ~= sr))
s2 %#ok
sr %#ok
error ('!') ;
end
err = norm (x-x3,1) ;
if (err > 1e-12)
x %#ok
x3 %#ok
uplo %#ok
err %#ok
error ('!')
end
maxerr = max (maxerr, err) ;
end
drawnow
end
fprintf ('maxerr = %g\n', maxerr) ;
|
