% Demo for colamd:  column 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:
% 
%	colamd		a replacement for colmmd, except that it doesn't
%			perform the column elimination tree post-ordering,
%			coletree, performed by colmmd.  Most similar in
%			function to the Matlab statement sparsfun ('colmmd', A).
%
%			Typical usage:  p = colamd (A) ;
%
%	colamdtree	colamd, followed by sparsfun ('coletree', A (:,p)).
%			Most similar in function to colmmd (A) itself.
%			For more information, do "type colmmd" in Matlab 5.2.
%
%			Typical usage:  p = colamdtree (A) ;
%
%	symamd		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).
%			Based on colamd.
%
%			Usage:  p = symamd (A) ;
%
%	symamdtree	symamd, followed by sparsfun ('symetree', A (p,p)).
%			Most similar in function to symmmd (A) itself.
%			Based on symamd.
%
%			Usage:  p = symamdtree (A) ;
%
% For a description of the methods used, see the colamd.c file.
%
% The following files form the colamd library:
%
%	colamd.README	a short abstract
%
%	colamd.c	the colamd computational kernel
%
%	colamd.h	include file for colamd.c, colamdmex.c, and symamdmex.c
%
%	colamd.m	a "dummy" m-file, with help information for colamd.
%			The actual ordering is done in a mexFunction.
%
%	colamd_demo.m	This file.
%
%	colamdmex.c	the colamd mexFunction.  Calls colamd.c.
%
%	colamdtree.m	the colamdtree m-file.  Calls colamd (A), then
%			sparsfun ('coletree', ...).
%
%	startup.m	running Matlab in the Colamd1.0/ directory automatically
%			runs this demo.
%
%	symamd.m	a "dummy" m-file, with help information for symamd.
%			The actual ordering is done in a mexFunction.
%
%	symamdmex.c	the symamd mexFunction.  Calls colamd.c
%
%	symamdtree.m	the symamdtree m-file.  Calls symamd (A), then
%			sparsfun ('symetree', ...).
%
% To compile the mexFunctions (the 2nd command is not a typo):
%
%	mex -O -output colamd colamd.c colamdmex.c
%	mex -O -output symamd colamd.c symamdmex.c
%
% (on Solaris 2.6, we use -xO4 -native, in the mexopts.sh file).
%
% Authors:
%
%	The authors of the code itself are Stefan I. Larimore and Timothy A.
%	Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was
%	developed in collaboration with John Gilbert, Xerox PARC, and Esmond
%	Ng, Oak Ridge National Laboratory.
%
% Date:
%
%	August 3, 1998.  Version 1.0.  Written under Matlab 5.2.0.3084, and
%	tested in Solaris 2.6.
%
% Acknowledgements:
%
%	This work was supported by the National Science Foundation, under
%	grants DMS-9504974 and DMS-9803599.
%
% Notice:
%
%	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 program 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 colamd" or "colamd"
%	(with no arguments) must cite the Authors.  Permission to modify the
%	code 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.
%
%	This software is provided free of charge.
%
% Availability:
%
%	http://www.cise.ufl.edu/~davis/colamd/
%	http://www.netlib.org/linalg/colamd/
%	http://www.mathworks.com/ as a user-contributed v5 package
%	http://www.netlib.org/scalapack/ colamd.c and colamd.h are part of the
%		SuperLU package

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

more on

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Colamd/symamd demo.') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help colamd_demo ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Colamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help colamd ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Colamdtree help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help colamdtree ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Symamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help symamd ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Symamdtree help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help symamdtree ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'colamd and symamd mexFunctions will now be compiled.') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, '\n\nHit any key to continue...\n') ;
pause
disp ('mex -O -output colamd colamd.c colamdmex.c')
mex -O -output colamd colamd.c colamdmex.c
disp ('mex -O -output symamd colamd.c symamdmex.c')
mex -O -output symamd colamd.c symamdmex.c

%-------------------------------------------------------------------------------
% 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) ;
b = (1:n)' ;

fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colamd:\n') ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	q = colamd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I = speye (n) ;
Q = I (:, q) ;
flops (0) ;
[L,U,P] = lu (A*Q) ;
x = Q * (U \ (L \ (P * b))) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q) and\n') ;
fprintf (1, 'then x = Q * (U \\ (L \\ (P * b))):           %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 lu (PAQ = LU), where Q is from colmmd:\n') ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	q = colmmd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I = speye (n) ;
Q = I (:, q) ;
flops (0) ;
[L,U,P] = lu (A*Q) ;
x = Q * (U \ (L \ (P * b))) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q) and\n') ;
fprintf (1, 'then x = Q * (U \\ (L \\ (P * b))):           %d\n', flops) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

fprintf (1, '\n\nTesting backslash - colmmd turned OFF:\n') ;
spparms ('autommd', 0) ;
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\nRestoring default spparms (this is recommended):\n') ;
spparms ('default') ;
spparms

%-------------------------------------------------------------------------------
% Large demo for colamd
%-------------------------------------------------------------------------------

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

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

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

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

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = colamd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = colmmd (A) ;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%-------------------------------------------------------------------------------
% Large demo for symamd
%-------------------------------------------------------------------------------

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

fprintf (1, 'Generating a random symmetric %d-by-%d sparse matrix.\n', n, 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\nSymamd run time:\n') ;
tic ;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	p = symamd (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

.