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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
|
%CCOLAMD_DEMO demo for ccolamd and csymamd
% minimum degree ordering algorithm.
%
% Example:
% ccolamd_demo
%
% See also ccolamd
% CCOLAMD, Copyright (c) 2005-2022, Univ. of Florida, All Rights Reserved.
% Authors: Timothy A. Davis, Sivasankaran Rajamanickam, and Stefan Larimore.
% SPDX-License-Identifier: BSD-3-clause
%-------------------------------------------------------------------------------
% Print the introduction, the help info, and compile the mexFunctions
%-------------------------------------------------------------------------------
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'ccolamd/csymamd demo.') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help ccolamd_demo ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'ccolamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help ccolamd ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'csymamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help csymamd ;
%-------------------------------------------------------------------------------
% Solving Ax=b
%-------------------------------------------------------------------------------
n = 100 ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Solving Ax=b for a small %d-by-%d random matrix:', n, n) ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, '\nNote: Random sparse matrices are AWFUL test cases.\n') ;
fprintf (1, 'They''re just easy to generate in a demo.\n') ;
% set up the system
rand ('state', 0) ;
randn ('state', 0) ;
spparms ('default') ;
A = sprandn (n, n, 2/n) + speye (n) ;
b = (1:n)' ;
clf ;
subplot (3,4,1)
spy (A)
title ('original matrix')
fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from ccolamd:\n') ;
q = ccolamd (A, 1) ;
I = speye (n) ;
Q = I (:, q) ;
[L,U,P] = lu (A*Q) ;
fl = luflops (L, U) ;
x = Q * (U \ (L \ (P * b))) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ;
fprintf (1, 'residual: %e\n', norm (A*x-b));
subplot (3,4,2) ;
spy (L|U) ;
title ('LU with ccolamd') ;
try
fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colamd:\n') ;
q = colamd (A) ;
I = speye (n) ;
Q = I (:, q) ;
[L,U,P] = lu (A*Q) ;
fl = luflops (L, U) ;
x = Q * (U \ (L \ (P * b))) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ;
fprintf (1, 'residual: %e\n', norm (A*x-b));
subplot (3,4,3) ;
spy (L|U) ;
title ('LU with colamd') ;
catch
fprintf (1, 'You have a very old version of MATLAB (no colamd) \n') ;
end
fprintf (1, '\n\nSolving via lu (PA = LU), without regard for sparsity:\n') ;
[L,U,P] = lu (A) ;
fl = luflops (L, U) ;
x = U \ (L \ (P * b)) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ;
fprintf (1, 'residual: %e\n', norm (A*x-b));
subplot (3,4,4) ;
spy (L|U) ;
title ('LU with no ordering') ;
%-------------------------------------------------------------------------------
% Large demo for ccolamd
%-------------------------------------------------------------------------------
% Since the analysis will be done on the Cholesky factorization of A'A,
% set the knob to tell ccolamd to order for Cholesky, not LU.
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Large demo for ccolamd (symbolic analysis only):') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
rand ('state', 0) ;
randn ('state', 0) ;
spparms ('default') ;
n = 1000 ;
fprintf (1, 'Generating a random %d-by-%d sparse matrix.\n', n, n) ;
A = sprandn (n, n, 2/n) + speye (n) ;
subplot (3,4,5)
spy (A)
title ('original matrix')
fprintf (1, '\n\nUnordered matrix:\n') ;
[lnz,h,parent,post,R] = symbfact (A, 'col') ;
fprintf (1, 'nz in Cholesky factors of A''A: %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A''A: %d\n', sum (lnz.^2)) ;
subplot (3,4,6) ;
spy (R) ;
title ('Cholesky with no ordering') ;
tic ;
p = ccolamd (A) ;
t = toc ;
[lnz,h,parent,post,R] = symbfact (A (:,p), 'col') ;
fprintf (1, '\n\nccolamd run time: %f\n', t) ;
fprintf (1, 'ccolamd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p): %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', sum (lnz.^2)) ;
subplot (3,4,7) ;
spy (R) ;
title ('Cholesky with ccolamd') ;
try
tic ;
p = colamd (A) ;
t = toc ;
[lnz,h,parent,post,R] = symbfact (A (:,p), 'col') ;
fprintf (1, '\n\ncolamd run time: %f\n', t) ;
fprintf (1, 'colamd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p): %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', sum (lnz.^2)) ;
subplot (3,4,8) ;
spy (R) ;
title ('Cholesky with colamd') ;
catch
fprintf (1, 'You have a very old version of MATLAB (no colamd) \n') ;
end
%-------------------------------------------------------------------------------
% Large demo for csymamd
%-------------------------------------------------------------------------------
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Large demo for csymamd (symbolic analysis only):') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Generating a random symmetric %d-by-%d sparse matrix.\n', n, n) ;
A = A+A' ;
subplot (3,4,9) ;
spy (A)
title ('original matrix')
fprintf (1, '\n\nUnordered matrix:\n') ;
[lnz,h,parent,post,R] = symbfact (A, 'sym') ;
fprintf (1, 'nz in Cholesky factors of A: %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A: %d\n', sum (lnz.^2)) ;
subplot (3,4,10) ;
spy (R) ;
title ('Cholesky with no ordering') ;
tic ;
p = csymamd (A) ;
t = toc ;
[lnz,h,parent,post,R] = symbfact (A (p,p), 'sym') ;
fprintf (1, '\n\ncsymamd run time: %f\n', t) ;
fprintf (1, 'csymamd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(p,p): %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;
subplot (3,4,11) ;
spy (R) ;
title ('Cholesky with csymamd') ;
try
tic ;
p = symamd (A) ;
t = toc ;
lnz = symbfact (A (p,p), 'sym') ;
fprintf (1, '\n\nsymamd run time: %f\n', t) ;
fprintf (1, 'symamd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(p,p): %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;
subplot (3,4,12) ;
spy (R) ;
title ('Cholesky with symamd') ;
catch
fprintf (1, 'You have a very old version of MATLAB (no symamd) \n') ;
end
drawnow
|
