% Demo for amdbar:  approximate minimum degree ordering algorithm.
% 
% The following m-files and mexFunctions provide alternative sparse matrix
% ordering methods for Matlab.  They are typically faster (sometimes much
% faster) and typically provide better orderings than their Matlab counterparts:
%
%	amdbar		a replacement for symmmd, except that it doesn't
%			perform the symmetric elimination tree post-ordering,
%			symetree, performed by symmmd.  Most similar in
%			function to the Matlab statement sparsfun ('symmmd', A).
%			Also similar to symamd (A).
%
% You need the following files:
%
%	amdbarmex.f	the Matlab interface to amdbar 
%
%	amdbar.f	the AMDBAR computational kernel
%
%	amdbar.m	so that "help amdbar" prints something useful
%
%	amdbar_demo.m	This file.
%
% To compile the mexFunction:
%
%	mex -O -output amdbar amdbarmex.f amdbar.f
%
% Authors:
%
%	Amdbar was written by Timothy A. Davis, Patrick Amestoy, Iain S. Duff,
%	and John K. Reid.  Timothy A. Davis (davis@cise.ufl.edu), University
%	of Florida, wrote the Matlab interface for amdbar.
%
% Date (of the Matlab interface for AMDBAR):
%
%	August 6, 1998.  Version 1.0.  
%
% Acknowledgements:
%
%	This work was supported by the National Science Foundation, under
%	grants DMS-9504974 and DMS-9803599.
%
% Notice (note the difference between amdbarmex.f and amdbar.f, below):
%
%	Copyright (c) 1998 by the University of Florida.  All Rights Reserved.
%
%	THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
%	EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
%
%	Permission is hereby granted to use or copy this file (ONLY) for any
%	purpose, provided the above notices are retained on all copies.
%	User documentation of any code that uses this code must cite the
%	Authors, the Copyright, and "Used by permission."  If this code is
%	accessible from within Matlab, then typing "help amdbar"
%	(with no arguments) must cite the Authors.  Permission to modify this
%	file (ONLY) and to distribute modified code is granted, provided the
%	above notices are retained, and a notice that the code was modified is
%	included with the above copyright notice.  You must also retain the
%	Availability information below, of the original version.
%
%	NOTE: This Matlab interface software is provided free of charge.
%	However, the computational kernel, amdbar.f, has stricter licensing
%	requirements.  Please see the licensing restrictions in amdbar.f,
%	specifically:
%
%	************************************************************************
%	* NOTICE:  "The AMD routines (AMDEXA, AMDBAR, AMDHAF, AMDHAT, AMDTRU,
%	* and AMDATR) may be used SOLELY for educational, research, and
%	* benchmarking purposes by non-profit organizations and the U.S.
%	* government.  Commercial and other organizations may make use of the
%	* AMD routines SOLELY for benchmarking purposes only.  The AMD
%	* routines may be modified by or on behalf of the User for such
%	* use but at no time shall the AMD routines or any such modified
%	* version of them become the property of the User.  The AMD routines
%	* are provided without warranty of any kind, either expressed or
%	* implied.  Neither the Authors nor their employers shall be liable
%	* for any direct or consequential loss or damage whatsoever arising
%	* out of the use or misuse of the AMD routines by the User.  The AMD
%	* routines must not be sold.  You may make copies of the AMD routines,
%	* but this NOTICE and the Copyright notice must appear in all copies.
%	* Any other use of the AMD routines requires written permission.
%	* Your use of the AMD routines is an implicit agreement to these
%	* conditions."
%	************************************************************************
%
%	You *must* abide by these restrictions for amdbar.f, even though this
%	file (amdbarmex.f, the Matlab C	interface for AMDBAR) has less stringent
%	restrictions.
%
% Availability:
%
%	http://www.cise.ufl.edu/~davis/amd/
%	http://www.netlib.org/linalg/amd/
%	For MC47BD, a faster version of AMDBAR:  Harwell
%	Subroutine Library, B 552, AEA Technology, Harwell, Didcot,
%	Oxon OX11 0RA, England.   telephone (44) 1235 434573,
%	fax (44) 1235 434340, email john.harding@aeat.co.uk, who will provide
%	details of price and conditions of use.

%-------------------------------------------------------------------------------
% Print the introduction, the help info, and compile the mexFunctions
%-------------------------------------------------------------------------------

more on

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'amdbar demo.') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help amdbar_demo ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'amdbar mexFunction will now be compiled.') ;
fprintf (1, '\n(ignore the warning, if any, about the variable "dext").') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, '\n\nHit any key to continue...\n') ;
pause
disp ('mex -O -output amdbar amdbarmex.f amdbar.f')
mex -O -output amdbar amdbarmex.f amdbar.f

%-------------------------------------------------------------------------------
% Solving Ax=b
%-------------------------------------------------------------------------------

n = 100 ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Solving Ax=b for a small %d-by-%d random matrix:', n, n) ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, '\nNOTE: Random sparse matrices are AWFUL test cases.\n') ;
fprintf (1, 'They''re just easy to generate in a demo.\n') ;

% set up the system

rand ('state', 0) ;
randn ('state', 0) ;
A = sprandn (n, n, 5/n) + speye (n) ;
A = A'*A ;
b = (1:n)' ;

fprintf (1, '\n\nSolving via chol (PAP'' = LL''), where P is from amdbar:\n') ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = amdbar (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

flops (0) ;
R = chol (A (p,p)) ;
x = R \ (R' \ b (p)) ;
x (p) = x ;
fprintf (1, '\nFlop count for R = chol (A (p,p)) and\n') ;
fprintf (1, 'then x (p) = R \\ (R'' \\ b (p)) :           %d\n', flops) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

fprintf (1, '\n\nTesting backslash - uses colmmd by default:\n') ;
spparms ('default') ;
spparms
flops (0) ;
x = A\b ;
fprintf (1, '\nFlop count for x=A\\b:                       %d\n', flops) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

fprintf (1, '\n\nSolving via chol (PAP'' = LL''), where P is from symmmd:\n') ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = symmmd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

flops (0) ;
R = chol (A (p,p)) ;
x = R \ (R' \ b (p)) ;
x (p) = x ;
fprintf (1, '\nFlop count for R = chol (A (p,p)) and\n') ;
fprintf (1, 'then x (p) = R \\ (R'' \\ b (p)) :           %d\n', flops) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

%-------------------------------------------------------------------------------
% Large demo for amdbar
%-------------------------------------------------------------------------------

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Demo for amdbar:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;

rand ('state', 0) ;
randn ('state', 0) ;
spparms ('default') ;
n = 10000 ;
fprintf (1, 'Generating a random symmetric %d-by-%d sparse matrix.\n', n, n) ;
A = sprandn (n, n, 5/n) + speye (n) ;
A = A + A' ;

fprintf (1, '\n\nUnordered matrix:\n') ;
lnz = symbfact (A, 'sym') ;
fprintf (1, 'nz in Cholesky factors of A:       %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A:      %d\n', sum (lnz.^2)) ;

fprintf (1, '\n\namdbar run time:\n') ;
tic ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = amdbar (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

toc ;
lnz = symbfact (A (p,p), 'sym') ;
fprintf (1, 'symamd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(p,p):  %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;

fprintf (1, '\n\nsymmmd run time:\n') ;
tic ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = symmmd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

toc ;
lnz = symbfact (A (p,p), 'sym') ;
fprintf (1, 'symmmd ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(p,p):  %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;

more off

.