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