% Most matrices occurring in practice
% have special structure.
% 1. Sparse Matrices
% MATLAB has facilities for handling sparse
% matrices, matrices with few nonzero elements.
a = sparse(4,6,3,10,10);
% we defined a 10-by-10 matrix with one
% nonzero element a(4,6) = 3
a

a =

   (4,6)        3

spy(a) % provides a picture
% Sierpinski gasket
p32 = pascal(32);
s32 = rem(p32,2); % remove the even entries
spy(s32) % patterns shows the Sierpinski gasket
32*32

ans =

        1024

p64 = pascal(64);
s64 = rem(p64,2); % remove the even entries
spy(s64) % patterns shows the Sierpinski gasket
%
%  It is interested to see the "hole" in the gasket,
%  due to the exponentional growth in the numbers,
%  MATLAB is no longer able to determine the parity correctly.
%
help sparfun
  Sparse matrices.
 
  Elementary sparse matrices.
    <a href="matlab:help speye">speye</a>       - Sparse identity matrix.
    <a href="matlab:help sprand">sprand</a>      - Sparse uniformly distributed random matrix.
    <a href="matlab:help sprandn">sprandn</a>     - Sparse normally distributed random matrix.
    <a href="matlab:help sprandsym">sprandsym</a>   - Sparse random symmetric matrix.
    <a href="matlab:help spdiags">spdiags</a>     - Sparse matrix formed from diagonals.
 
  Full to sparse conversion.
    <a href="matlab:help sparse">sparse</a>      - Create sparse matrix.
    <a href="matlab:help full">full</a>        - Convert sparse matrix to full matrix.
    <a href="matlab:help find">find</a>        - Find indices of nonzero elements.
    <a href="matlab:help spconvert">spconvert</a>   - Import from sparse matrix external format.
 
  Working with sparse matrices.
    <a href="matlab:help nnz">nnz</a>         - Number of nonzero matrix elements.
    <a href="matlab:help nonzeros">nonzeros</a>    - Nonzero matrix elements.
    <a href="matlab:help nzmax">nzmax</a>       - Amount of storage allocated for nonzero matrix elements.
    <a href="matlab:help spones">spones</a>      - Replace nonzero sparse matrix elements with ones.
    <a href="matlab:help spalloc">spalloc</a>     - Allocate space for sparse matrix.
    <a href="matlab:help issparse">issparse</a>    - True for sparse matrix.
    <a href="matlab:help spfun">spfun</a>       - Apply function to nonzero matrix elements.
    <a href="matlab:help spy">spy</a>         - Visualize sparsity pattern.
 
  Reordering algorithms.
    <a href="matlab:help colamd">colamd</a>      - Column approximate minimum degree permutation.
    <a href="matlab:help symamd">symamd</a>      - Symmetric approximate minimum degree permutation.
    <a href="matlab:help colmmd">colmmd</a>      - Column minimum degree permutation.
    <a href="matlab:help symmmd">symmmd</a>      - Symmetric minimum degree permutation.
    <a href="matlab:help symrcm">symrcm</a>      - Symmetric reverse Cuthill-McKee permutation.
    <a href="matlab:help colperm">colperm</a>     - Column permutation.
    <a href="matlab:help randperm">randperm</a>    - Random permutation.
    <a href="matlab:help dmperm">dmperm</a>      - Dulmage-Mendelsohn permutation.
 
  Linear algebra.
    <a href="matlab:help eigs">eigs</a>        - A few eigenvalues, using ARPACK.
    <a href="matlab:help svds">svds</a>        - A few singular values, using eigs.
    <a href="matlab:help luinc">luinc</a>       - Incomplete LU factorization.
    <a href="matlab:help cholinc">cholinc</a>     - Incomplete Cholesky factorization.
    <a href="matlab:help normest">normest</a>     - Estimate the matrix 2-norm.
    <a href="matlab:help condest">condest</a>     - 1-norm condition number estimate.
    <a href="matlab:help sprank">sprank</a>      - Structural rank.
 
  Linear Equations (iterative methods).
    <a href="matlab:help pcg">pcg</a>         - Preconditioned Conjugate Gradients Method.
    <a href="matlab:help bicg">bicg</a>        - BiConjugate Gradients Method.
    <a href="matlab:help bicgstab">bicgstab</a>    - BiConjugate Gradients Stabilized Method.
    <a href="matlab:help cgs">cgs</a>         - Conjugate Gradients Squared Method.
    <a href="matlab:help gmres">gmres</a>       - Generalized Minimum Residual Method.
    <a href="matlab:help lsqr">lsqr</a>        - Conjugate Gradients on the Normal Equations.
    <a href="matlab:help minres">minres</a>      - Minimum Residual Method.
    <a href="matlab:help qmr">qmr</a>         - Quasi-Minimal Residual Method.
    <a href="matlab:help symmlq">symmlq</a>      - Symmetric LQ Method.
 
  Operations on graphs (trees).
    <a href="matlab:help treelayout">treelayout</a>  - Lay out tree or forest.
    <a href="matlab:help treeplot">treeplot</a>    - Plot picture of tree.
    <a href="matlab:help etree">etree</a>       - Elimination tree.
    <a href="matlab:help etreeplot">etreeplot</a>   - Plot elimination tree.
    <a href="matlab:help gplot">gplot</a>       - Plot graph, as in "graph theory".
 
  Miscellaneous.
    <a href="matlab:help symbfact">symbfact</a>    - Symbolic factorization analysis.
    <a href="matlab:help spparms">spparms</a>     - Set parameters for sparse matrix routines.
    <a href="matlab:help spaugment">spaugment</a>   - Form least squares augmented system.



% Let us convert a sparse matrix into dense format
s8 = rem(pascal(8),2);
size(8,8);
size(s8)

ans =

     8     8

s8a = sparse(s8)

s8a =

   (1,1)        1
   (2,1)        1
   (3,1)        1
   (4,1)        1
   (5,1)        1
   (6,1)        1
   (7,1)        1
   (8,1)        1
   (1,2)        1
   (3,2)        1
   (5,2)        1
   (7,2)        1
   (1,3)        1
   (2,3)        1
   (5,3)        1
   (6,3)        1
   (1,4)        1
   (5,4)        1
   (1,5)        1
   (2,5)        1
   (3,5)        1
   (4,5)        1
   (1,6)        1
   (3,6)        1
   (1,7)        1
   (2,7)        1
   (1,8)        1

[i,j,s] = find(s8)

i =

     1
     2
     3
     4
     5
     6
     7
     8
     1
     3
     5
     7
     1
     2
     5
     6
     1
     5
     1
     2
     3
     4
     1
     3
     1
     2
     1


j =

     1
     1
     1
     1
     1
     1
     1
     1
     2
     2
     2
     2
     3
     3
     3
     3
     4
     4
     5
     5
     5
     5
     6
     6
     7
     7
     8


s =

     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1
     1

% in the triplet (i,j,s) we find the row
% indices in i, column indices in j, and
% the values at (i,j) in s
% These triplets are useful to define sparse
% matrices:
sp8b = sparse(i,j,s);
spy(sp8b);
% to define a 50-by-50 tridiagonal matrix,
% with 2 on the diagonal, and ones on the
% lower and upper diagonal
i = [ 1:50  2:50  1:49 ];
j = [ 1:50  1:49  2:50 ];
s = [ 2*ones(1,50)  ones(1,49)  ones(1,49) ];
d3 = sparse(i,j,s);
spy(d3)
d3(13,20:30)

ans =

   All zero sparse: 1-by-11

d3(25,20:30)

ans =

   (1,5)        1
   (1,6)        2
   (1,7)        1

% 8.2 is another application
% 8.3 Multi-Dimensional matrices
r = -2:0.1:2;
[x,y,z] = ndgrid(r,r,r);
t = y.*exp(-x.^2 - y.^2 - z.^2);
slice(y,x,z,t,[],[-1.2  .8],-2);
diary off