%% Sparse Matrices
% This M-file considers full and sparse matrices
% It looks at their construction, manipulation and visualization.
%%
% A sparse matrix is a matrix where the majority of the entries are zero,
% in Matlab a special type of matrix is used which "squeezes" out the zero
% entries (it saves memory).
%%
% This is a random matrix:
N=5;
M=rand(N);
M2=M*M; % This is its matrix square.
Mp2=M.*M; % This is its point-wise square.
%%
% This is an identity matrix and a sparse identity matrix:
myI=eye(N)
myspI=speye(N)
%%
% I can check if a matrix is sparse like this:
issparse(myI)
issparse(myspI)
%%
% I can force a matrix to be sparse:
mysamp=3*myI; mysamp(2,3)=-2; 
mysamp=sparse(mysamp);
%%
% This looks dumb but can be quite effective especially if we want to use
% the sparse matrix many times later.
%%
% This is an example of a sparse matrix build in action.
% It builds a second order differentiation matrix.
%%
% Define a spatial grid:
xmin = 0;
xmax = 50;
N = 2^7;
x = linspace(xmin,xmax,N+1);
dx=x(2)-x(1); dx2=dx*dx;
%%
% Now build a vector of ones,
e=ones(N+1,1);
%%
% Next use spdiags to drop the correct vectors into the correct diagonals.
Dxx = spdiags([e -2*e e], -1:1, N+1, N+1);
% finally scale
Dxx=(1/dx2)*Dxx;
%%
% This is how you visualize a matrix structure. Filled entries are in blue
% while empty ones are white.
figure(1)
clf
spy(Dxx);
%%
% Sparsity can be very helpful when carrying out matrix solves since it
% offers immense speed up especially for very large matrices.
Dxxf=full(Dxx); % This converts Dxx to a full matrix.
myrhs=sin(8*pi*x/(xmax-xmin))';
tic,mysols=Dxx\myrhs; toc
tic,mysolsf=Dxxf\myrhs; toc
max(abs(mysols-mysolsf))

