Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/MATLAB/umfpack_demo.m @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 function umfpack_demo
2 % UMFPACK DEMO
3 %
4 % A demo of UMFPACK for MATLAB.
5 %
6 % See also umfpack, umfpack_make, umfpack_details, umfpack_report,
7 % and umfpack_simple.
8
9 % UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.
10 % 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 MATLAB:\n') ;
45 xm = A\b ;
46 x = xm ;
47
48 fprintf ('Difference between UMFPACK and MATLAB 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 spy (A)
60 title ('The matrix A') ;
61
62 subplot (2,3,2)
63 [P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;
64 treeplot (Fr (1:end-1,2)') ;
65 title ('Supernodal column elimination tree') ;
66
67 subplot (2,3,3)
68 spy (P1 * A * Q1)
69 title ('A, with initial row and column order') ;
70
71 subplot (2,3,4)
72 fprintf ('\n--------------------------------------------------------------\n') ;
73 fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ;
74 [L, U, P, Q, R] = umfpack (A) ;
75 spy (P*A*Q)
76 title ('A, with final row/column order') ;
77
78 fprintf ('\nP * (R\\A) * Q - L*U should be zero:\n') ;
79 fprintf ('norm (P*(R\\A)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ...
80 norm (P * (R\A) * Q - L*U, 1), lu_normest (P * (R\A) * Q, L, U)) ;
81
82 fprintf ('\nSolution to Ax=b via UMFPACK factorization:\n') ;
83 fprintf ('x = Q * (U \\ (L \\ (P * (R \\ b))))\n') ;
84 xu = Q * (U \ (L \ (P * (R \ b)))) ;
85 x = xu ;
86
87 fprintf ('\nUMFPACK flop count: %d\n', luflop (L, U)) ;
88
89 subplot (2,3,5)
90 spy (spones (L) + spones (U))
91 title ('UMFPACK LU factors') ;
92
93 subplot (2,3,6)
94 fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\n') ;
95 fprintf ('If you are using a version of MATLAB prior to V6.0, then the\n') ;
96 fprintf ('following statement (q = colamd (A)) may fail. Either download\n');
97 fprintf ('colamd from http://www.cise.ufl.edu/research/sparse, upgrade to\n') ;
98 fprintf ('MATLAB V6.0 or later, or replace the statement with\n') ;
99 fprintf ('q = colmmd (A) ;\n') ;
100 try
101 q = colamd (A) ;
102 catch
103 fprintf ('\n *** colamd not found, using colmmd instead *** \n') ;
104 q = colmmd (A) ;
105 end
106 [L, U, P] = lu (A (:,q)) ;
107 spy (spones (L) + spones (U))
108 title ('MATLAB LU factors') ;
109
110 fprintf ('\nSolution to Ax=b via MATLAB factorization:\n') ;
111 fprintf ('x = U \\ (L \\ (P * b)) ; x (q) = x ;\n') ;
112 xm = U \ (L \ (P * b)) ;
113 xm (q) = xm ;
114
115 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
116 norm (xu - xm, Inf)) ;
117
118 fprintf ('\nMATLAB LU flop count: %d\n', luflop (L, U)) ;
119
120 %-------------------------------------------------------------------------------
121 % solve A'x=b
122 %-------------------------------------------------------------------------------
123
124 fprintf ('\n--------------------------------------------------------------\n') ;
125 fprintf ('Solve A''x=b:\n') ;
126
127 fprintf ('Solving A''x=b via UMFPACK:\n') ;
128 [xu, info] = umfpack (b', '/', A, control) ;
129 xu = xu' ;
130
131 fprintf ('Solving A''x=b via MATLAB:\n') ;
132 xm = (b'/A)' ;
133 x = xm ;
134
135 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
136 norm (xu - xm, Inf)) ;
137
138 %-------------------------------------------------------------------------------
139 % factor A' and then solve Ax=b using the factors of A'
140 %-------------------------------------------------------------------------------
141
142 fprintf ('\n--------------------------------------------------------------\n') ;
143 fprintf ('Compute C = A'', and compute the LU factorization of C.\n') ;
144 fprintf ('Factorizing A'' can sometimes be better than factorizing A itself\n');
145 fprintf ('(less work and memory usage). Solve C''x=b; the solution is the\n') ;
146 fprintf ('same as the solution to Ax=b for the original A.\n');
147
148 C = A' ;
149
150 % factorize C (P,Q) = L*U
151 [L, U, P, Q, R, info] = umfpack (C, control) ;
152
153 fprintf ('\nP * (R\\C) * Q - L*U should be zero:\n') ;
154 fprintf ('norm (P*(R\\C)*Q - L*U, 1) = %g (exact) %g (estimated)\n', ...
155 norm (P * (R\C) * Q - L*U, 1), lu_normest (P * (R\C) * Q, L, U)) ;
156
157 fprintf ('\nSolution to Ax=b via UMFPACK, using the factors of C:\n') ;
158 fprintf ('x = R \\ (P'' * (L'' \\ (U'' \\ (Q'' * b)))) ;\n') ;
159 xu = R \ (P' * (L' \ (U' \ (Q' * b)))) ;
160 x = xu ;
161
162 fprintf ('Solution to Ax=b via MATLAB:\n') ;
163 xm = A\b ;
164 x = xm ;
165
166 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
167 norm (xu - xm, Inf)) ;
168
169 %-------------------------------------------------------------------------------
170 % solve Ax=B
171 %-------------------------------------------------------------------------------
172
173 fprintf ('\n--------------------------------------------------------------\n') ;
174 fprintf ('\nSolve AX=B, where B is n-by-10, and sparse\n') ;
175 B = sprandn (n, 10, 0.05) ;
176 XU = umfpack_solve (A, '\', B, control) ;
177 XM = A\B ;
178
179 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
180 norm (XU - XM, Inf)) ;
181
182 fprintf ('\n--------------------------------------------------------------\n') ;
183 fprintf ('\nSolve AX=B, where B is n-by-10, and sparse, using umfpack_btf\n') ;
184 XU = umfpack_btf (A, B, control) ;
185
186 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
187 norm (XU - XM, Inf)) ;
188
189 fprintf ('\n--------------------------------------------------------------\n') ;
190 fprintf ('\nSolve A''X=B, where B is n-by-10, and sparse\n') ;
191 XU = umfpack_solve (B', '/', A, control) ;
192 XM = B'/A ;
193
194 fprintf ('Difference between UMFPACK and MATLAB solution: %g\n', ...
195 norm (XU - XM, Inf)) ;
196
197 %-------------------------------------------------------------------------------
198 % compute the determinant
199 %-------------------------------------------------------------------------------
200
201 fprintf ('\n--------------------------------------------------------------\n') ;
202 fprintf ('det(A): %g UMFPACK determinant: %g\n', det (A), umfpack (A, 'det')) ;