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
|
function test21
%TEST21 test cs_updown, chol_updown2
%
% Example:
% test21
% See also: testall
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
clear functions
rand ('state', 0) ;
randn ('state', 0) ;
clf
for trials = 1:100
if (trials <= 1)
n = trials ;
else
n = 1+fix (100 * rand (1)) ;
end
fprintf ('n: %d\n', n) ;
d = 0.1 * rand (1) ;
A = sprandn (n,n,d) ;
if (~ispc)
if (rand ( ) > .5)
A = A + 1i * sprand (A) ;
end
end
A = A+A' + 100 * speye (n) ;
try
p = amd (A) ;
catch
p = symamd (A) ;
end
A = sparse (A (p,p)) ;
try
L = chol (A)' ;
catch
continue ;
end
parent = etree (A) ;
subplot (1,3,1) ;
spy (A) ;
if (n > 0)
subplot (1,3,2) ;
treeplot (parent) ;
end
subplot (1,3,3) ;
spy (L) ;
drawnow
for trials2 = 1:10
k = 1+fix (n * rand (1)) ;
if (k <= 0 | k > n) %#ok
k = 1 ;
end
w = sprandn (L (:,k)) ;
Anew = A + w*w' ;
Lnew = cs_updown (L, w, parent) ;
err6 = norm (Lnew*Lnew' - Anew, 1) ;
Lnew = cs_updown (L, w, parent, '+') ;
err7 = norm (Lnew*Lnew' - Anew, 1) ;
[Lnew, wnew] = chol_updown2 (L, 1, w) ;
err2 = norm (Lnew*Lnew' - Anew, 1) ;
err10 = norm (wnew - (L\w)) ;
L3 = chol_updown2 (L, +1, w) ;
err9 = norm (L3*L3' - Anew, 1) ;
[L2, wnew] = chol_updown2 (Lnew, -1, w) ;
err3 = norm (L2*L2' - A, 1) ;
err11 = norm (wnew - (Lnew\w)) ;
L2 = cs_updown (Lnew, w, parent, '-') ;
err5 = norm (L2*L2' - A, 1) ;
L2 = chol_updown2 (Lnew, -1, w) ;
err8 = norm (L2*L2' - A, 1) ;
err = max ([err2 err3 err5 err6 err7 err9 err8 err10 err11]) ;
fprintf (' k %3d %6.2e\n', k, err) ;
if (err > 1e-11)
err2 %#ok
err3 %#ok
err5 %#ok
err6 %#ok
err7 %#ok
err8 %#ok
err9 %#ok
err10 %#ok
err11 %#ok
pause
end
end
% pause
end
|
