--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/OCTAVE/umfpack_demo.m	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,191 @@
+function umfpack_demo
+% UMFPACK DEMO
+%
+% A demo of UMFPACK for OCTAVE.
+%
+% See also umfpack, umfpack_make, umfpack_details, umfpack_report,
+% and umfpack_simple.
+
+% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 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 OCTAVE:\n') ;
+xm = A\b ;
+x = xm ;
+
+fprintf ('Difference between UMFPACK and OCTAVE solution: %g\n', ...
+    norm (xu - xm, Inf)) ;
+
+%-------------------------------------------------------------------------------
+% spy the results
+%-------------------------------------------------------------------------------
+
+figure (1) ;
+clf
+
+subplot (2,3,1)
+title ('The matrix A') ;
+spy (A)
+
+subplot (2,3,2)
+[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;
+title ('Supernodal column elimination tree') ;
+%% Disable for now !!
+%% treeplot (Fr (1:end-1,2)') ;
+
+subplot (2,3,3)
+title ('A, with initial row and column order') ;
+spy (P1 * A * Q1)
+
+subplot (2,3,4)
+fprintf ('\n--------------------------------------------------------------\n') ;
+fprintf ('\nFactorizing [L, U, P, Q, R] = umfpack (A)\n') ;
+[L, U, P, Q, R] = umfpack (A) ;
+title ('A, with final row/column order') ;
+spy (P*A*Q)
+
+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)
+title ('UMFPACK LU factors') ;
+spy (spones (L) + spones (U))
+
+subplot (2,3,6)
+fprintf ('\nFactorizing [L, U, P] = lu (A (:, q))\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)) ;
+title ('OCTAVE LU factors') ;
+spy (spones (L) + spones (U))
+
+fprintf ('\nSolution to Ax=b via OCTAVE 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 OCTAVE solution: %g\n', ...
+    norm (xu - xm, Inf)) ;
+
+fprintf ('\nOCTAVE 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 OCTAVE:\n') ;
+xm = (b'/A)' ;
+x = xm ;
+
+fprintf ('Difference between UMFPACK and OCTAVE 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 OCTAVE:\n') ;
+xm = A\b ;
+x = xm ;
+
+fprintf ('Difference between UMFPACK and OCTAVE 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 OCTAVE 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 OCTAVE 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 OCTAVE solution: %g\n', ...
+    norm (XU - XM, Inf)) ;