lxsd / fvn

Blame view

fvn_sparse/AMD/MATLAB/amd_demo.m 2.38 KB
edit raw normal view history
422234dc3   daniau   git-svn-id: https... 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
 
  function amd_demo
  %AMD_DEMO a demo of amd2, using the can_24 matrix
  %
  % A demo of AMD for MATLAB.
  %
  % Example:
  %   amd_demo
  %
  % See also: amd, amd2, amd_make
  
  % Copyright 1994-2007, Tim Davis, University of Florida,
  % Patrick R. Amestoy, and Iain S. Duff. 
  
  % This orders the same matrix as the ANSI C demo, amd_demo.c.  It includes an
  % additional analysis of the matrix via MATLAB's symbfact routine.
  
  % First, print the help information for AMD
  help amd2
  
  % Get the Harwell/Boeing can_24 matrix.
  
  load can_24
  A = spconvert (can_24) ;
  
  n = size (A,1) ;
  
  figure (1)
  clf
  hold off
  subplot (2,2,1) ;
  spy (A)
  title ('HB/can24 matrix') ;
  
  % order the matrix.  Note that the Info argument is optional.
  fprintf ('
  If the next step fails, then you have
  ') ;
  fprintf ('not yet compiled the AMD mexFunction.
  ') ;
  [p, Info] = amd2 (A) ;		%#ok
  
  % order again, but this time print some statistics
  [p, Info] = amd2 (A, [10 1 1]) ;
  
  fprintf ('Permutation vector:
  ') ;
  fprintf (' %2d', p) ;
  fprintf ('
  
  ') ;
  
  subplot (2,2,2) ;
  spy (A (p,p)) ;
  title ('Permuted matrix') ;
  
  % The amd_demo.c program stops here.
  
  fprintf ('Analyze A(p,p) with MATLAB''s symbfact routine:
  ') ;
  [cn, height, parent, post, R] = symbfact (A (p,p)) ;
  
  subplot (2,2,3) ;
  spy (R') ; 
  title ('Cholesky factor, L') ;
  
  subplot (2,2,4) ;
  treeplot (parent) ;
  title ('elimination tree') ;
  
  % results from symbfact
  lnz = sum (cn) ;                % number of nonzeros in L, incl. diagonal
  cn = cn - 1 ;                   % get the count of off-diagonal entries
  fl = n + sum (cn.^2 + 2*cn) ;   % flop count for chol (A (p,p)
  fprintf ('number of nonzeros in L (including diagonal):      %d
  ', lnz) ;
  fprintf ('floating point operation count for chol (A (p,p)): %d
  ', fl) ;
  
  % approximations from amd:
  lnz2 = n + Info (10) ;
  fl2 = n + Info (11) + 2 * Info (12) ;
  fprintf ('
  Results from AMD''s approximate analysis:
  ') ;
  fprintf ('number of nonzeros in L (including diagonal):      %d
  ', lnz2) ;
  fprintf ('floating point operation count for chol (A (p,p)): %d
  
  ', fl2) ;
  
  if (lnz2 ~= lnz | fl ~= fl2)						    %#ok
      fprintf ('Note that the nonzero and flop counts from AMD are slight
  ') ;
      fprintf ('upper bounds.  This is due to the approximate minimum degree
  ');
      fprintf ('method used, in conjunction with "mass elimination".
  ') ;
      fprintf ('See the discussion about mass elimination in amd.h and
  ') ;
      fprintf ('amd_2.c for more details.
  ') ;
  end