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
|
function test7
%TEST7 test cs_lu
%
% Example:
% test7
% See also: testall
% CSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
index = ssget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;
f = f (1:100) ;
clf
for i = f
Prob = ssget (i) ;
disp (Prob) ;
A = Prob.A ;
if (~isreal (A))
continue
end
[m n] = size (A) ;
if (m ~= n)
continue
end
[L,U,P] = lu (A) ;
udiag = full (diag (U)) ;
umin = min (abs (udiag)) ;
fprintf ('umin %g\n', umin) ;
if (umin > 1e-14)
[L2,U2,p] = cs_lu (A) ;
subplot (3,4,1) ; spy (A) ;
subplot (3,4,2) ; spy (A(p,:)) ;
subplot (3,4,3) ; spy (L2) ;
subplot (3,4,4) ; spy (U2) ;
err1 = norm (L*U-P*A,1) ;
err2 = norm (L2*U2-A(p,:),1) ;
fprintf ('err %g %g\n', err1, err2) ;
end
q = colamd (A) ;
[L,U,P] = lu (A (:,q)) ;
udiag = full (diag (U)) ;
umin = min (abs (udiag)) ;
fprintf ('umin %g with q\n', umin) ;
if (umin > 1e-14)
[L2,U2,p,q2] = cs_lu (A) ;
subplot (3,4,5) ; spy (A) ;
subplot (3,4,6) ; spy (A(p,q2)) ;
subplot (3,4,7) ; spy (L2) ;
subplot (3,4,8) ; spy (U2) ;
err1 = norm (L*U-P*A(:,q),1) ;
err2 = norm (L2*U2-A(p,q2),1) ;
fprintf ('err %g %g\n', err1, err2) ;
end
try
q = amd (A) ;
catch
q = symamd (A) ;
end
tol = 0.01 ;
[L,U,P] = lu (A (q,q), tol) ;
udiag = full (diag (U)) ;
umin = min (abs (udiag)) ;
fprintf ('umin %g with amd q\n', umin) ;
if (umin > 1e-14)
[L2,U2,p,q2] = cs_lu (A,tol) ;
subplot (3,4,9) ; spy (A) ;
subplot (3,4,10) ; spy (A(p,q2)) ;
subplot (3,4,11) ; spy (L2) ;
subplot (3,4,12) ; spy (U2) ;
err1 = norm (L*U-P*A(q,q),1) ;
err2 = norm (L2*U2-A(p,q2),1) ;
lbig = full (max (max (abs (L2)))) ;
fprintf ('err %g %g lbig %g\n', err1, err2, lbig) ;
if (lbig > 1/tol)
error ('L!') ;
end
end
drawnow
% pause
end
|
