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