home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / SPARSITY.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.9 KB  |  146 lines
  1. %SPARSITY Demonstrate effect of sparsity orderings.
  2.  
  3. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  4.  
  5. clf
  6. clc
  7. echo on
  8.  
  9. % This demonstration shows the effect of different orderings
  10. % in a sparse matrix on the time and storage requirements for
  11. % sparse matrix operations.
  12. % We begin with a Harwell-Boeing sparse test matrix, west0479.
  13. % This matrix describes the connections in a model of a
  14. % diffraction column in a chemical plant.
  15.  
  16. load('west0479.mat')
  17. A = west0479;
  18.  
  19. % A "spy" plot shows the locations of the nonzero elements.
  20.  
  21. spy(A)
  22. title('spy(A)')
  23.  
  24. % Press any key to continue after pauses.
  25. pause
  26.  
  27. clc
  28.  
  29. % Here are spy plots of A, A', A*A' and A'*A.
  30.  
  31. subplot(2,2,1), spy(A),    title('A'),     drawnow
  32. subplot(2,2,2), spy(A'),   title('A'''),   drawnow
  33. subplot(2,2,3), spy(A*A'), title('A*A'''), drawnow
  34. subplot(2,2,4), spy(A'*A); title('A''*A'), drawnow
  35.  
  36. % The spy plots of A*A' and A'*A show "fill-in", the creation of
  37. % nonzero elements in positions that are zero in the original matrix.
  38. % Storage requirements for sparse matrices are proportional to
  39. % the number of nonzero elements.
  40.  
  41. pause
  42.  
  43. clf
  44. clc
  45.  
  46. % We continue with A*A', which is symmetric and positive definite.
  47.  
  48. S = A*A';
  49.  
  50. subplot(2,2,1), spy(S), title('S')
  51.  
  52. pause
  53.  
  54. clc
  55.  
  56. % The computation of the Cholesky factorization creates fill-in.
  57. % Let's time the factorization and count the number of nonzeros
  58. % in the triangular factor.
  59.  
  60. tic
  61. R = chol(S);
  62. toc
  63. nz = nnz(R)
  64.  
  65. subplot(2,2,2), spy(R), title('chol(S)')
  66.  
  67. pause
  68.  
  69. clc
  70.  
  71. % The reverse Cuthill-McKee algorithm renumbers the rows and columns
  72. % in an attempt to reduce the bandwidth or profile of the matrix.
  73.  
  74. p = symrcm(S);
  75. subplot(2,2,3), spy(S(p,p)), title('S(p,p)')
  76.  
  77. pause
  78.  
  79. clc
  80.  
  81. % The fill-in is confined to the band, so the Cholesky factorization
  82. % of the reordered matrix takes less time and less storage.
  83.  
  84. tic
  85. R = chol(S(p,p));
  86. toc
  87. nz = nnz(R)
  88.  
  89. subplot(2,2,4), spy(R), title('chol(S(p,p))')
  90.  
  91. pause
  92.  
  93. clf
  94. clc
  95.  
  96. % The column count ordering moves rows and columns with higher
  97. % nonzero count towards the end of the matrix.
  98.  
  99. q = colperm(S);
  100. subplot(2,2,1), spy(S(q,q)), title('S(q,q)')
  101.  
  102. pause
  103.  
  104. clc
  105.  
  106. % For this example, the column count ordering happens to reduce
  107. % the time and storage for Cholesky factorization, but this
  108. % behavior cannot be expected in general.
  109.  
  110. tic
  111. R = chol(S(q,q));
  112. toc
  113. nz = nnz(R)
  114.  
  115. subplot(2,2,2), spy(R), title('chol(S(q,q))')
  116.  
  117. pause
  118.  
  119. clc
  120.  
  121. % Minimum degree reordering is a powerful, graph theoretic algorithm
  122. % that produces large blocks of zeros in the matrix.
  123.  
  124. r = symmmd(S);
  125. subplot(2,2,3), spy(S(r,r)), title('S(r,r)')
  126.  
  127. pause
  128.  
  129. clc
  130.  
  131. % These blocks are preserved during the Cholesky factorization.
  132. % As a result, minimum degree requires the least time and storage.
  133.  
  134. tic
  135. R = chol(S(r,r));
  136. toc
  137. nz = nnz(R)
  138.  
  139. subplot(2,2,4), spy(R), title('chol(S(r,r))')
  140.  
  141. pause
  142.  
  143. % End
  144. echo off
  145. disp('>>')
  146.