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
|
function umfpack_demo
% UMFPACK DEMO
%
% A demo of UMFPACK for MATLAB.
%
% See also umfpack, umfpack_make, umfpack_details, umfpack_report,
% and umfpack_simple.
% UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.
% All Rights Reserved. Type umfpack_details for License.
%-------------------------------------------------------------------------------
% get default control parameters
%-------------------------------------------------------------------------------
control = umfpack ;
fprintf ('\nEnter the printing level for UMFPACK''s output statistics:\n') ;
fprintf ('0: none, 1: errors only, 2: statistics, 4: print some of outputs\n') ;
c = input ('5: print all output [default is 1]: ') ;
if (isempty (c))
c = 1 ;
end
control (1) = c ;
%-------------------------------------------------------------------------------
% solve a simple system
%-------------------------------------------------------------------------------
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('Factor and solve a small system, Ax=b, using default parameters\n') ;
if (control (1) > 1)
fprintf ('(except for verbose printing enabled)\n') ;
end
load west0067
A = Problem.A ;
n = size (A, 1) ;
b = rand (n, 1) ;
fprintf ('Solving Ax=b via UMFPACK:\n') ;
[xu, info] = umfpack (A, '\', b, control) ;
x = xu ;
fprintf ('Solving Ax=b via MATLAB:\n') ;
xm = A\b ;
x = xm ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (xu - xm, Inf)) ;
%-------------------------------------------------------------------------------
% spy the results
%-------------------------------------------------------------------------------
figure (1)
clf
subplot (2,3,1)
spy (A)
title ('The matrix A') ;
subplot (2,3,2)
[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;
treeplot (Fr (1:end-1,2)') ;
title ('Supernodal column elimination tree') ;
subplot (2,3,3)
spy (P1 * A * Q1)
title ('A, with initial row and column order') ;
subplot (2,3,4)
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ;
[L, U, P, Q, R] = umfpack (A) ;
spy (P*A*Q)
title ('A, with final row/column order') ;
fprintf ('\nP * (R\\A) * Q - L*U should be zero:\n') ;
fprintf ('norm (P*(R\\A)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ...
norm (P * (R\A) * Q - L*U, 1), lu_normest (P * (R\A) * Q, L, U)) ;
fprintf ('\nSolution to Ax=b via UMFPACK factorization:\n') ;
fprintf ('x = Q * (U \\ (L \\ (P * (R \\ b))))\n') ;
xu = Q * (U \ (L \ (P * (R \ b)))) ;
x = xu ;
fprintf ('\nUMFPACK flop count: %d\n', luflop (L, U)) ;
subplot (2,3,5)
spy (spones (L) + spones (U))
title ('UMFPACK LU factors') ;
subplot (2,3,6)
fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\n') ;
fprintf ('If you are using a version of MATLAB prior to V6.0, then the\n') ;
fprintf ('following statement (q = colamd (A)) may fail. Either download\n');
fprintf ('colamd from http://www.cise.ufl.edu/research/sparse, upgrade to\n') ;
fprintf ('MATLAB V6.0 or later, or replace the statement with\n') ;
fprintf ('q = colmmd (A) ;\n') ;
try
q = colamd (A) ;
catch
fprintf ('\n *** colamd not found, using colmmd instead *** \n') ;
q = colmmd (A) ;
end
[L, U, P] = lu (A (:,q)) ;
spy (spones (L) + spones (U))
title ('MATLAB LU factors') ;
fprintf ('\nSolution to Ax=b via MATLAB factorization:\n') ;
fprintf ('x = U \\ (L \\ (P * b)) ; x (q) = x ;\n') ;
xm = U \ (L \ (P * b)) ;
xm (q) = xm ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (xu - xm, Inf)) ;
fprintf ('\nMATLAB LU flop count: %d\n', luflop (L, U)) ;
%-------------------------------------------------------------------------------
% solve A'x=b
%-------------------------------------------------------------------------------
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('Solve A''x=b:\n') ;
fprintf ('Solving A''x=b via UMFPACK:\n') ;
[xu, info] = umfpack (b', '/', A, control) ;
xu = xu' ;
fprintf ('Solving A''x=b via MATLAB:\n') ;
xm = (b'/A)' ;
x = xm ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (xu - xm, Inf)) ;
%-------------------------------------------------------------------------------
% factor A' and then solve Ax=b using the factors of A'
%-------------------------------------------------------------------------------
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('Compute C = A'', and compute the LU factorization of C.\n') ;
fprintf ('Factorizing A'' can sometimes be better than factorizing A itself\n');
fprintf ('(less work and memory usage). Solve C''x=b; the solution is the\n') ;
fprintf ('same as the solution to Ax=b for the original A.\n');
C = A' ;
% factorize C (P,Q) = L*U
[L, U, P, Q, R, info] = umfpack (C, control) ;
fprintf ('\nP * (R\\C) * Q - L*U should be zero:\n') ;
fprintf ('norm (P*(R\\C)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ...
norm (P * (R\C) * Q - L*U, 1), lu_normest (P * (R\C) * Q, L, U)) ;
fprintf ('\nSolution to Ax=b via UMFPACK, using the factors of C:\n') ;
fprintf ('x = R \\ (P'' * (L'' \\ (U'' \\ (Q'' * b)))) ;\n') ;
xu = R \ (P' * (L' \ (U' \ (Q' * b)))) ;
x = xu ;
fprintf ('Solution to Ax=b via MATLAB:\n') ;
xm = A\b ;
x = xm ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (xu - xm, Inf)) ;
%-------------------------------------------------------------------------------
% solve Ax=B
%-------------------------------------------------------------------------------
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('\nSolve AX=B, where B is n-by-10, and sparse\n') ;
B = sprandn (n, 10, 0.05) ;
XU = umfpack_solve (A, '\', B, control) ;
XM = A\B ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (XU - XM, Inf)) ;
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('\nSolve AX=B, where B is n-by-10, and sparse, using umfpack_btf\n') ;
XU = umfpack_btf (A, B, control) ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (XU - XM, Inf)) ;
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('\nSolve A''X=B, where B is n-by-10, and sparse\n') ;
XU = umfpack_solve (B', '/', A, control) ;
XM = B'/A ;
fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
norm (XU - XM, Inf)) ;
%-------------------------------------------------------------------------------
% compute the determinant
%-------------------------------------------------------------------------------
fprintf ('\n--------------------------------------------------------------\n') ;
fprintf ('det(A): %g UMFPACK determinant: %g\n', det (A), umfpack (A, 'det')) ;
|
