← Back to branch summary

~logan/ubuntu/trusty/suitesparse/4.2.1-3ubuntu1

« back to all changes in this revision

Viewing changes to LDL/MATLAB/ldldemo.m

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-05-29 09:36:29 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070529093629-zowquo0b7slkk6nc
Tags: 3.0.0-2
http://bugs.debian.org/426349
http://bugs.debian.org/426354
* suitesparse builds properly twice in a row
* Bug fix: "suitesparse - FTBFS: Broken build depens: libgfortran1-dev",
  thanks to Bastian Blank (Closes: #426349).
* Bug fix: "suitesparse_3.0.0-1: FTBFS: build-depends on
  libgfortran1-dev", thanks to Steve Langasek (Closes: #426354).

Show diffs side-by-side

added added

removed removed

Lines of Context:
LDL/MATLAB/ldldemo.m
 
1
function ldldemo
 
2
%LDLDEMO demo program for LDL
 
3
%
 
4
% Example:
 
5
%   ldldemo
 
6
%
 
7
% See also ldlsparse.
 
8
 
 
9
% Copyright 2006-2007 by Timothy A. Davis, Univ. of Florida
 
10
 
 
11
% compile the LDLSPARSE and LDLSYMBOL mexFunctions
 
12
help ldlsparse
 
13
 
 
14
fprintf ('\nTesting ldlsparse and ldlsymbol:\n') ;
 
15
 
 
16
% create a small random symmetric positive definite sparse matrix
 
17
n = 100 ;
 
18
d = 0.03 ;
 
19
rand ('state', 0) ;
 
20
randn ('state', 0) ;
 
21
A = sprandn (n, n, d) ;
 
22
A = speye (n) + A*A' ;
 
23
b = randn (n, 1) ;
 
24
 
 
25
figure (1)
 
26
clf
 
27
subplot (2,2,1) ;
 
28
spy (A) ;
 
29
title ('original matrix') ;
 
30
 
 
31
% permute for sparsity
 
32
p = symamd (A) ;
 
33
C = A (p,p) ;
 
34
 
 
35
subplot (2,2,2) ;
 
36
spy (C) ;
 
37
title ('permuted matrix') ;
 
38
drawnow
 
39
 
 
40
% factorize, without using ldlsparse's internal permutation
 
41
[L, D, Parent, fl] = ldlsparse (C) ;
 
42
L = L + speye (n) ;
 
43
err = norm (L*D*L' - C, 1) ;
 
44
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;
 
45
 
 
46
% solve Ax=b
 
47
x = L' \ (D \ (L \ (b (p)))) ;
 
48
x (p) = x ;
 
49
resid = norm (A*x-b) ;
 
50
fprintf ('residual %g for ldlsparse, flops %10.1f\n', resid, fl) ;
 
51
 
 
52
% solve Ax=b with one call to ldlsparse
 
53
x = ldlsparse (C, [ ], b (p)) ;
 
54
x (p) = x ;
 
55
resid = norm (A*x-b) ;
 
56
fprintf ('residual %g for ldlsparse solve\n', resid) ;
 
57
 
 
58
subplot (2,2,3) ;
 
59
spy (L + D + L') ;
 
60
title ('L+D+L''') ;
 
61
 
 
62
subplot (2,2,4) ;
 
63
treeplot (Parent)
 
64
title ('elimination tree') ;
 
65
 
 
66
% try ldlrow (this will be slow)
 
67
[L, D] = ldlrow (C) ;
 
68
x = L' \ (D \ (L \ (b (p)))) ;
 
69
x (p) = x ;
 
70
resid = norm (A*x-b) ;
 
71
fprintf ('residual %g for ldlrow.m\n', resid) ;
 
72
 
 
73
% factorize, using ldlsparse's internal permutation
 
74
[L, D, Parent, fl] = ldlsparse (A, p) ;
 
75
L = L + speye (n) ;
 
76
err = norm (L*D*L' - C, 1) ;
 
77
fprintf ('norm (LDL''-PAP'') = %g\n', err) ;
 
78
 
 
79
% solve Ax=b
 
80
x = L' \ (D \ (L \ (b (p)))) ;
 
81
x (p) = x ;
 
82
resid = norm (A*x-b) ;
 
83
fprintf ('residual %g for ldlsparse, flops %10.1f\n', resid, fl) ;
 
84
 
 
85
% solve Ax=b with one call to ldlsparse
 
86
x = ldlsparse (A, p, b) ;
 
87
resid = norm (A*x-b) ;
 
88
fprintf ('residual %g for ldlsparse solve\n\n', resid) ;
 
89
 
 
90
% compare ldlsymbol and symbfact
 
91
[Lnz, Parent, fl] = ldlsymbol (A) ;
 
92
fprintf ('Original matrix: nz in L: %5d  flop count: %g\n', sum (Lnz), fl) ;
 
93
 
 
94
Lnz2 = symbfact (A) - 1 ;
 
95
Parent2 = etree (A) ;
 
96
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
 
97
if (any (Lnz ~= Lnz2))
 
98
    error ('Lnz mismatch') ;
 
99
end
 
100
if (any (Parent ~= Parent2))
 
101
    error ('Parent mismatch') ;
 
102
end
 
103
if (fl ~= fl2)
 
104
    error ('fl mismatch') ;
 
105
end
 
106
 
 
107
[Lnz, Parent, fl] = ldlsymbol (A, p) ;
 
108
fprintf ('Permuted matrix: nz in L: %5d  flop count: %g\n', sum (Lnz), fl) ;
 
109
 
 
110
Lnz2 = symbfact (A (p,p)) - 1 ;
 
111
Parent2 = etree (A (p,p)) ;
 
112
fl2 = sum (Lnz2 .* (Lnz2 + 2)) ;
 
113
if (any (Lnz ~= Lnz2))
 
114
    error ('Lnz mismatch') ;
 
115
end
 
116
if (any (Parent ~= Parent2))
 
117
    error ('Parent mismatch') ;
 
118
end
 
119
if (fl ~= fl2)
 
120
    error ('fl mismatch') ;
 
121
end
 
122
 
 
123
 
 
124
fprintf ('\nldldemo: all tests passed\n') ;