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
|
%
% ldldemo.m: demo program for the LDL and LDLSYMBOL mexFunctions
%
% LDL Version 1.2, Copyright (c) 2005 by Timothy A Davis,
% University of Florida. All Rights Reserved. See README for the License.
% compile the LDL and LDLSYMBOL mexFunctions
clear
help ldl
fprintf ('Compiling ldl and ldlsymbol:\n') ;
mex -inline ldl.c ldlmex.c
mex -inline -output ldlsymbol ldl.c ldlsymbolmex.c
fprintf ('\nTesting ldl and ldlsymbol:\n') ;
% create a small random symmetric positive definite sparse matrix
n = 100 ;
d = 0.03 ;
rand ('state', 0) ;
randn ('state', 0) ;
A = sprandn (n, n, d) ;
A = speye (n) + A*A' ;
b = randn (n, 1) ;
figure (1)
clf
subplot (2,2,1) ;
spy (A) ;
title ('original matrix') ;
% permute for sparsity
p = symamd (A) ;
C = A (p,p) ;
subplot (2,2,2) ;
spy (C) ;
title ('permuted matrix') ;
drawnow
% factorize, without using ldl's internal permutation
[L, D, Parent, fl] = ldl (C) ;
L = L + speye (n) ;
err = norm (L*D*L' - C, 1) ;
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;
% solve Ax=b
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldl, flops %10.1f\n', resid, fl) ;
% solve Ax=b with one call to ldl
x = ldl (C, [ ], b (p)) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldl solve\n', resid) ;
subplot (2,2,3) ;
spy (L + D + L') ;
title ('L+D+L''') ;
subplot (2,2,4) ;
treeplot (Parent)
title ('elimination tree') ;
% try ldlrow (this will be slow)
[L, D] = ldlrow (C) ;
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldlrow.m\n', resid) ;
% factorize, using ldl's internal permutation
[L, D, Parent, fl] = ldl (A, p) ;
L = L + speye (n) ;
err = norm (L*D*L' - C, 1) ;
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;
% solve Ax=b
x = L' \ (D \ (L \ (b (p)))) ;
x (p) = x ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldl, flops %10.1f\n', resid, fl) ;
% solve Ax=b with one call to ldl
x = ldl (A, p, b) ;
resid = norm (A*x-b) ;
fprintf ('residual %g for ldl solve\n\n', resid) ;
% compare ldlsymbol and symbfact
[Lnz, Parent, fl] = ldlsymbol (A) ;
fprintf ('Original matrix: nz in L: %5d flop count: %g\n', sum (Lnz), fl) ;
Lnz2 = symbfact (A) - 1 ;
Parent2 = etree (A) ;
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
if (any (Lnz ~= Lnz2)) error ('Lnz mismatch') ; end
if (any (Parent ~= Parent2)) error ('Parent mismatch') ; end
if (fl ~= fl2) error ('fl mismatch') ; end
[Lnz, Parent, fl] = ldlsymbol (A, p) ;
fprintf ('Permuted matrix: nz in L: %5d flop count: %g\n', sum (Lnz), fl) ;
Lnz2 = symbfact (A (p,p)) - 1 ;
Parent2 = etree (A (p,p)) ;
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
if (any (Lnz ~= Lnz2)) error ('Lnz mismatch') ; end
if (any (Parent ~= Parent2)) error ('Parent mismatch') ; end
if (fl ~= fl2) error ('fl mismatch') ; end
fprintf ('\nldldemo: all tests passed\n') ;
|
