|
5164
|
1 function umfpack_demo |
|
|
2 % UMFPACK DEMO |
|
|
3 % |
|
|
4 % A demo of UMFPACK for OCTAVE. |
|
|
5 % |
|
|
6 % See also umfpack, umfpack_make, umfpack_details, umfpack_report, |
|
|
7 % and umfpack_simple. |
|
|
8 |
|
|
9 % UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
|
10 % Davis. All Rights Reserved. Type umfpack_details for License. |
|
|
11 |
|
|
12 %------------------------------------------------------------------------------- |
|
|
13 % get default control parameters |
|
|
14 %------------------------------------------------------------------------------- |
|
|
15 |
|
|
16 control = umfpack ; |
|
|
17 fprintf ('\nEnter the printing level for UMFPACK''s output statistics:\n') ; |
|
|
18 fprintf ('0: none, 1: errors only, 2: statistics, 4: print some of outputs\n') ; |
|
|
19 c = input ('5: print all output [default is 1]: ') ; |
|
|
20 if (isempty (c)) |
|
|
21 c = 1 ; |
|
|
22 end |
|
|
23 control (1) = c ; |
|
|
24 |
|
|
25 %------------------------------------------------------------------------------- |
|
|
26 % solve a simple system |
|
|
27 %------------------------------------------------------------------------------- |
|
|
28 |
|
|
29 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
30 fprintf ('Factor and solve a small system, Ax=b, using default parameters\n') ; |
|
|
31 if (control (1) > 1) |
|
|
32 fprintf ('(except for verbose printing enabled)\n') ; |
|
|
33 end |
|
|
34 |
|
|
35 load west0067 ; |
|
|
36 A = Problem.A ; |
|
|
37 n = size (A, 1) ; |
|
|
38 b = rand (n, 1) ; |
|
|
39 |
|
|
40 fprintf ('Solving Ax=b via UMFPACK:\n') ; |
|
|
41 [xu, info] = umfpack (A, '\\', b, control) ; |
|
|
42 x = xu ; |
|
|
43 |
|
|
44 fprintf ('Solving Ax=b via OCTAVE:\n') ; |
|
|
45 xm = A\b ; |
|
|
46 x = xm ; |
|
|
47 |
|
|
48 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
49 norm (xu - xm, Inf)) ; |
|
|
50 |
|
|
51 %------------------------------------------------------------------------------- |
|
|
52 % spy the results |
|
|
53 %------------------------------------------------------------------------------- |
|
|
54 |
|
|
55 figure (1) ; |
|
|
56 clf |
|
|
57 |
|
|
58 subplot (2,3,1) |
|
|
59 title ('The matrix A') ; |
|
|
60 spy (A) |
|
|
61 |
|
|
62 subplot (2,3,2) |
|
|
63 [P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ; |
|
|
64 title ('Supernodal column elimination tree') ; |
|
|
65 %% Disable for now !! |
|
|
66 %% treeplot (Fr (1:end-1,2)') ; |
|
|
67 |
|
|
68 subplot (2,3,3) |
|
|
69 title ('A, with initial row and column order') ; |
|
|
70 spy (P1 * A * Q1) |
|
|
71 |
|
|
72 subplot (2,3,4) |
|
|
73 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
74 fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ; |
|
|
75 [L, U, P, Q, R] = umfpack (A) ; |
|
|
76 title ('A, with final row/column order') ; |
|
|
77 spy (P*A*Q) |
|
|
78 |
|
|
79 fprintf ('\nP * (R\\A) * Q - L*U should be zero:\n') ; |
|
|
80 fprintf ('norm (P*(R\\A)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... |
|
|
81 norm (P * (R\A) * Q - L*U, 1), lu_normest (P * (R\A) * Q, L, U)) ; |
|
|
82 |
|
|
83 fprintf ('\nSolution to Ax=b via UMFPACK factorization:\n') ; |
|
|
84 fprintf ('x = Q * (U \\ (L \\ (P * (R \\ b))))\n') ; |
|
|
85 xu = Q * (U \ (L \ (P * (R \ b)))) ; |
|
|
86 x = xu ; |
|
|
87 |
|
|
88 fprintf ('\nUMFPACK flop count: %d\n', luflop (L, U)) ; |
|
|
89 |
|
|
90 subplot (2,3,5) |
|
|
91 title ('UMFPACK LU factors') ; |
|
|
92 spy (spones (L) + spones (U)) |
|
|
93 |
|
|
94 subplot (2,3,6) |
|
|
95 fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\n') ; |
|
|
96 try |
|
|
97 q = colamd (A) ; |
|
|
98 catch |
|
|
99 fprintf ('\n *** colamd not found, using colmmd instead *** \n') ; |
|
|
100 q = colmmd (A) ; |
|
|
101 end |
|
|
102 [L, U, P] = lu (A (:,q)) ; |
|
|
103 title ('OCTAVE LU factors') ; |
|
|
104 spy (spones (L) + spones (U)) |
|
|
105 |
|
|
106 fprintf ('\nSolution to Ax=b via OCTAVE factorization:\n') ; |
|
|
107 fprintf ('x = U \\ (L \\ (P * b)) ; x (q) = x ;\n') ; |
|
|
108 xm = U \ (L \ (P * b)) ; |
|
|
109 xm (q) = xm ; |
|
|
110 |
|
|
111 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
112 norm (xu - xm, Inf)) ; |
|
|
113 |
|
|
114 fprintf ('\nOCTAVE LU flop count: %d\n', luflop (L, U)) ; |
|
|
115 |
|
|
116 %------------------------------------------------------------------------------- |
|
|
117 % solve A'x=b |
|
|
118 %------------------------------------------------------------------------------- |
|
|
119 |
|
|
120 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
121 fprintf ('Solve A''x=b:\n') ; |
|
|
122 |
|
|
123 fprintf ('Solving A''x=b via UMFPACK:\n') ; |
|
|
124 [xu, info] = umfpack (b', '/', A, control) ; |
|
|
125 xu = xu' ; |
|
|
126 |
|
|
127 fprintf ('Solving A''x=b via OCTAVE:\n') ; |
|
|
128 xm = (b'/A)' ; |
|
|
129 x = xm ; |
|
|
130 |
|
|
131 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
132 norm (xu - xm, Inf)) ; |
|
|
133 |
|
|
134 %------------------------------------------------------------------------------- |
|
|
135 % factor A' and then solve Ax=b using the factors of A' |
|
|
136 %------------------------------------------------------------------------------- |
|
|
137 |
|
|
138 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
139 fprintf ('Compute C = A'', and compute the LU factorization of C.\n') ; |
|
|
140 fprintf ('Factorizing A'' can sometimes be better than factorizing A itself\n'); |
|
|
141 fprintf ('(less work and memory usage). Solve C''x=b; the solution is the\n') ; |
|
|
142 fprintf ('same as the solution to Ax=b for the original A.\n'); |
|
|
143 |
|
|
144 C = A' ; |
|
|
145 |
|
|
146 % factorize C (P,Q) = L*U |
|
|
147 [L, U, P, Q, R, info] = umfpack (C, control) ; |
|
|
148 |
|
|
149 fprintf ('\nP * (R\\C) * Q - L*U should be zero:\n') ; |
|
|
150 fprintf ('norm (P*(R\\C)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ... |
|
|
151 norm (P * (R\C) * Q - L*U, 1), lu_normest (P * (R\C) * Q, L, U)) ; |
|
|
152 |
|
|
153 fprintf ('\nSolution to Ax=b via UMFPACK, using the factors of C:\n') ; |
|
|
154 fprintf ('x = R \\ (P'' * (L'' \\ (U'' \\ (Q'' * b)))) ;\n') ; |
|
|
155 xu = R \ (P' * (L' \ (U' \ (Q' * b)))) ; |
|
|
156 x = xu ; |
|
|
157 |
|
|
158 fprintf ('Solution to Ax=b via OCTAVE:\n') ; |
|
|
159 xm = A\b ; |
|
|
160 x = xm ; |
|
|
161 |
|
|
162 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
163 norm (xu - xm, Inf)) ; |
|
|
164 |
|
|
165 %------------------------------------------------------------------------------- |
|
|
166 % solve Ax=B |
|
|
167 %------------------------------------------------------------------------------- |
|
|
168 |
|
|
169 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
170 fprintf ('\nSolve AX=B, where B is n-by-10, and sparse\n') ; |
|
|
171 B = sprandn (n, 10, 0.05) ; |
|
|
172 XU = umfpack_solve (A, '\\', B, control) ; |
|
|
173 XM = A\B ; |
|
|
174 |
|
|
175 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
176 norm (XU - XM, Inf)) ; |
|
|
177 |
|
|
178 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
179 fprintf ('\nSolve AX=B, where B is n-by-10, and sparse, using umfpack_btf\n') ; |
|
|
180 XU = umfpack_btf (A, B, control) ; |
|
|
181 |
|
|
182 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
183 norm (XU - XM, Inf)) ; |
|
|
184 |
|
|
185 fprintf ('\n--------------------------------------------------------------\n') ; |
|
|
186 fprintf ('\nSolve A''X=B, where B is n-by-10, and sparse\n') ; |
|
|
187 XU = umfpack_solve (B', '/', A, control) ; |
|
|
188 XM = B'/A ; |
|
|
189 |
|
|
190 fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ... |
|
|
191 norm (XU - XM, Inf)) ; |
