% 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