%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
% (1) matrix * vector
%     MatVecCo
%     MatVecRo
% (2) matrix * matrix
%     MatMatDot
%     MatMatOuter
%     MatMatSax
%     MatMatVec
% (3) Matrix structure
%     MakeBlock
%     TridTril
%     ShowSparse
%     Circulant1
%     Circulant2
% (4) Diagnostics  
%     matBench
%     matflops
%     Prod_bound
%     spies
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
function  y = MatVecCO(A,x)
% y = MatVexCO(A,x)
% Computes the matrix-vector product y = A*x (via saxpys) where
% A is an mXn matrix and x is a column vector.
[m,n] = size(A);
y = zeros(m,1);
for j = 1:n
   y = y + A(:,j)*x(j);
end
%
% /$\       / --$----- \      -  
% |$|       | --$----- |      -  
% |$|   ___ | --$----- |      $  
% |$|   \   | --$----- |      -  
% |$| =  >  | --$----- |  *   -  
% |$|   /__ | --$----- |      -  
% |$|       | --$----- |      -  
% |$|       | --$----- |      -  
% \$/       \ --$----- /  
%-----------------------------------------------------------------------
function  y = MatVecRO(A,x)
% y = MatVexRO(A,x)
% Computes the matrix-vector product y = A*x (via saxpys) where
% A is an mXn matrix and x is a column vector.
[m,n] = size(A);
y = zeros(m,1);
for i = 1:m
   y(i) = A(i,:)*x;
end
%       / -------- \     /$\
%       | -------- |     |$|
%   $ = | $$$$$$$$ |     |$|
%       | -------- |  *  |$|
%       | -------- |     |$|
%       | -------- |     |$|
%       \ -------- /     \$
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
function C = MatMatDot(A,B)
% C = MatMatDot(A,B)
% computes matrix*matrix product C = A*B (via dot products) where
% A is an mXp matrix, B a pXn matrix.
[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j = 1:n
   % Compute j-th column of C
   for i = 1:m
      C(i,j) = A(i,:)*B(:,j);
   end
end
%     3  6             8          3  6
%  /--|---\   / -------- \     /--$---\
%  |--|---|   | -------- |     |--$---|
%  |--|---|   | -------- |     |--$---|
% 4|--$---| = | $$$$$$$$ |4    |--$---|
%  |--|---|   | -------- |  *  |--$---|
%  |--|---|   | -------- |     |--$---|
% 7\--|---/  7\ -------- /     |--$---|
%                             8\--$---/

%-----------------------------------------------------------------------
function C = MatMatOuter(A,B)
% C = MatMatOuter(A,B)
% computes matrix*matrix product C = A*B (via outer products) where
% A is an mXp matrix, B a pXn matrix.
[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for k = 1:p
   % Add-in k-th outer product
   C = C + A(:,k)*B(k,:);
end
%
%
%
%        6                   8             6
%  /------\         / ----$--- \     /------\
%  |------|   _8_   | ----$--- |     |------|
%  |------|   \     | ----$--- |     |------|
%  |------| =  >    | ----$--- |     |------|
%  |------|   /__   | ----$--- |  *  |$$$$$$|k=5
%  |------|   k=1   | ----$--- |     |------|
% 7\------/        7\ ----$--- /     |------|
%                         k=5       8\------/
%-----------------------------------------------------------------------
function C = MatMatSax(A,B)
% C = MatMatSax(A,B)
% computes matrix*matrix product C = A*B (via saxpys) where
% A is an mXp matrix, B a pXn matrix.
[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j = 1:n
   % Compute j-th column of C
   for k = 1:p
      C(:,j) = C(:,j) + A(:,k)*B(k,j);
   end
end
%        6                  8             6
%  /--$---\        / ----$--- \     /--.---\
%  |--$---|   _8_  | ----$--- |     |--.---|
%  |--$---|   \    | ----$--- |     |--.---|
%  |--$---| =  >   | ----$--- |     |--.---|
%  |--$---|   /__  | ----$--- |  * k|--$---|
%  |--$---|   k=1  | ----$--- |     |--.---|
% 7\--$---/       7\ ----$--- /     |--.---|
%     j=3                k=4       8\--.---/   
%                                      j=3
%     C              C              A             B
% here, a column of C is built by adding the contributions from
% all the columns of A, each multiplied by the corresponding element in B
%-----------------------------------------------------------------------
function C = MatMatVec(A,B)
% C = MatMatVec(A,B)
% computes matrix*matrix product C = A*B (via matrix-vector products) where
% A is an mXp matrix, B a pXn matrix.
[m,p] = size(A);
[p,n] = size(B);
C = zeros(m,n);
for j = 1:n
   % Compute j-th column of C
   C(:,j) = A*B(:,j);
%  this addition is unnecessary, since C contains zeros!
%  C(:,j) = C(:,j) + A*B(:,j);
end
%
%
%        6               8             6
%  /--$---\     / -------- \     /--$---\
%  |--$---|     | -------- |     |--$---|
%  |--$---|     | -------- |     |--$---|
%  |--$---| =   | -------- |     |--$---|
%  |--$---|     | -------- |  *  |--$---|
%  |--$---|     | -------- |     |--$---|
% 7\--$---/    7\ -------- /     |--$---|
%     j=3                       8\--$---/
%                                   j=
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
function A = MakeBlock(A_scalar,p)
% A = MakeBlock(A_scalar,p)
% A_scalar is an nXn matrix and p divides n. A is an (n/p)X(n/p)
% cell array that represents A_scalar as a block matrix with pXp blocks.
[n,n] = size(A_scalar);
m = n/p;
A = cell(m,m);
for i = 1:m
   for j = 1:m
      A{i,j} = A_scalar(1+(i-1)*p:i*p,1+(j-1)*p:j*p);
   end
end
%-----------------------------------------------------------------------
% function C = TridTril(A,B)
% A: tridiagonal matrix
% B: Lower triangular matrix
% n = input('enter size of problem:');
disp(sprintf('    n      f1      f2        f3      f4        t1        t2
 t3       t4'))
disp(sprintf('------------------------------------------------------------------
------------'))
for n = 8:8:80
v1 = rand(n-1,1);
v2 = rand(n,1);
v3 = rand(n-1,1);
A = diag(v1,-1)+diag(v2)+diag(v3,1);
B=tril(rand(n,n));
flops(0);tic; C1 = A*B ; f1 = flops;
t1=toc;
flops(0); tic;C2 = sparse(A)*sparse(B) ; f2 = flops;
t2=toc;
flops(0);tic;
C3=zeros(n,n);
   C3(1,:)=                  A(1,1)*B(1,:)+A(1,  2)*B(  2,:);
for j = 2:n-1
   C3(j,:)=A(j,j-1)*B(j-1,:)+A(j,j)*B(j,:)+A(j,j+1)*B(j+1,:);  
   end
   C3(n,:)=A(n,n-1)*B(n-1,:)+A(n,n)*B(n,:);
f3 = flops;t3=toc;
E3 = C3-C1;
error3 = norm(E3,inf);
flops(0);tic;
C4=zeros(n,n);
   C4(1,1:  2)=                      A(1,1)*B(1,1:  2)+A(1,  2)*B(  2,1:  2);
for j = 2:n-1
   C4(j,1:j+1)=A(j,j-1)*B(j-1,1:j+1)+A(j,j)*B(j,1:j+1)+A(j,j+1)*B(j+1,1:j+1);
end
   C4(n,1:n  )=A(n,n-1)*B(n-1,1:n  )+A(n,n)*B(n,1:n  );
f4 = flops;t4=toc;
E4 = C4-C1;
error4 = norm(E4,inf);
disp(sprintf(' %4.0f %7.0f   %5.0f   %7.0f   %5.0f  %8.3f  %8.3f  %8.3f  %8.3f',
n,f1,f2,f3,f4,t1,t2,t3,t4))
end  
%-----------------------------------------------------------------------
% script showSparse
%-----------------------------------------------------------------------
   i = 0;
for L = 100:100:1000
     i  = i+1;
   x(i) = L;
   A   = diag(2*ones(L,1)) - diag(ones(L-1,1),-1) - diag(ones(L-1,1),1);
   S_A = sparse(A);
   flops(0); y =   A*rand(L,1); f1(i) = flops;
   flops(0); y = S_A*rand(L,1); f2(i) = flops;
end
   semilogy(x,f1,x,f2)
   legend('Flops with full A','Flops with sparse A',0)
%-----------------------------------------------------------------------
% function C = Circulant1(a) : Element-wise
%  C  = / a1 a2 ... an   \
%       | an a1 ... an-1 |
%       .  ...  ...  ... .
%       \ a2 a3 ... a1   /
  function C = Circulant1(a)
% C = Circulant1(a) is a circulant matrix with first row equal to a
n = length(a);
C = zeros(n,n);
for i = 1:n
   for j = 1:n
      C(i,j) = a(rem(n-i+j,n)+1);
   end
end
%-----------------------------------------------------------------------
% function C = Circulant2(a): Vectorized, with array combinations
%  C  = / a1 a2 ... an   \
%       | an a1 ... an-1 |
%       .  ...  ...  ... .
%       \ a2 a3 ... a1   /
	function C=Circulant2(a)
% C = Circulant2(a) is a circulant matrix with first row equal to a
n = length(a);
C = zeros(n,n);
%caution: the following line will not work unless a is a row vector
C(1,:) = a;
for i = 2:n
   C(i,:) = [ C(i-1,n) C(i-1,1:n-1) ];
end
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
% Script File: matBench
%-----------------------------------------------------------------------
disp('    n         Dot         Saxpy       MatVec        Outer        Direct ')
disp('------------------------------------------------------------------------')
for L = [100,200,400]
   A = randn(L,L);
   B = randn(L,L);
   tic; C = MatMatDot(A,B); f1 = toc;
   tic; C = MatMatSax(A,B); f2 = toc;
   tic; C = MatMatVec(A,B); f3 = toc;
   tic; C = MatMatOuter(A,B); f4 = toc;
   tic; C = A*B; f5 = toc;
disp(sprintf('%6.0f %12.8f %12.8f %12.8f %12.8f %12.8f',L,f1,f2,f3,f4,f5))
end
%-----------------------------------------------------------------------
% script: matflops
%-----------------------------------------------------------------------
disp('    n         Dot         Saxpy       MatVec        Outer        Direct ')
disp('------------------------------------------------------------------------')
for L = [100,200,400]
   A = randn(L,L);
   B = randn(L,L);
   flops(0); C = MatMatDot(A,B); f1 = flops;
   flops(0); C = MatMatSax(A,B); f2 = flops;
   flops(0); C = MatMatVec(A,B); f3 = flops;
   flops(0); C = MatMatOuter(A,B); f4 = flops;
   flops(0); C = A*B; f5 = flops;
disp(sprintf('%6.0f %12.0f %12.0f %12.0f %12.0f %12.0f',L,f1,f2,f3,f4,f5))
end
%-----------------------------------------------------------------------
% Script File: ProdBound
% Examines the error in 3-digit matrix multiplication

clc
eps3 = .005          % 3-digit machine precision
nRepeat = 10;        % number of trials per n-value
disp('   n   1-norm factor    ')
disp('------------------------')
for n = 2:10
   s = 0;
   for r = 1:nRepeat
      A  = randn(n,n);
      B  = randn(n,n);
      A1 = A+eps3*(ones(n,n)-2*rand(n,n));
      B1 = B+eps3*(ones(n,n)-2*rand(n,n));
format short e
      C = A*B;
format long e
      E = C-A1*B1;
      s = s+norm(E,1)/(eps3*norm(A,1)*norm(B,1));
   end
   disp(sprintf('%4.0f    %8.3f   ',n,s/nRepeat))
en
%-----------------------------------------------------------------------
% script file: spies % a command for viewing matrices
   A=gallery('wathen',8,8);
subplot(1,2,1), spy(A), subplot(1,2,2), spy(chol(A))
figure
r=symrcm(A);
subplot(1,2,1), spy(A(r,r)), subplot(1,2,2), spy(chol(A(r,r)))
figure
m=symmmd(A);
subplot(1,2,1), spy(A(m,m)), subplot(1,2,2), spy(chol(A(m,m)))
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------

