%  <  Friday Feb.8 2002  >
%-----------------------------------------------------------------------
%------------------------
% function A = Hilb(n)
% just like hilb(n)
function A = Hilb(n)
A = zeros(n,n);
for i = 1:n
  %for j = 1:n
   for j = i:n
      A(i,j) = 1/(i+j-1);
      A(j,i) = A(i,j);
   end
end
% alternatively, using vectors:
function A = hilbV(n)
a=1./(1:(2*n-1));
for i = 1:n
   A(i,:) = a(i:(i+n-1));
end
%------------------------
%-----------------------------------------------------------------------
% script combinatorial
P = zeros(n,n);
P(:,1) = ones(n,1);
for i = 2:n
   for j = 2:i
      P(i,j) = P(i-1,j-1)+P(i-1,j);
   end
end
% lternatively, we can compute the array by using the fact that the entries
% in row i appear as coefficients in the expansion of (1+x)^(i-1):
% (1+x)^k = P(k,0)    + P(k,   1) x + P(k,   2) x^2 + ... + P(k,   k  ) x^k  or
%         = P{k+1}(1) + P{k+1}(2) x + P{k+1}(3) x^2 + ... + P{k+1}(k+1) x^k
% i.e. (1+x)^k --> P{k+1} = [a0 a1 a2 a3 ... ak], a vector of length (k+1)
% Then, since
% x * (a0 + a1*x+a2*x^2+...+ak*x^k) = (0+a0*x+a1*x^2+a2*x^3+...+ak*x^(k+1))
% x [a0 a1 a2 ... ak 0] = [0 a0 a1 a2 ... ak], a vector of length (k+2)
% so that 
%  since  (1+x)^(k+1) = [b0 b1 b2 ... bk b(k+1)] (length k+2)
% equals  (1+x)*(1+x)^k = (1+x)^k + x (1+x)^k =
%[a0 a1 a2 ... ak 0 ]+[0 a0 a1 a2 ... ak]=[a0 a0+a1 a1+a2 ... a(k-1)+ak ak]
% script combinatorialV
n=20;
P{1} = [1];
for i = 2:(n+1)   %  what happens if we use i=2:n+1 ?
  P{i} = [0 P{i-1}] + [P{i-1} 0];
end
%-----------------------------------------------------------------------
% 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               
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
%  <  Monday Feb.11 2002  >
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
Banded Matrices
%-----------------------------------------------------------------------
% script tridiag
A=[0 1 2 3;-1 0 1 2;-2 -1 0 1;-3 -2 -1 0;-4 -3 -2 -1];
v=diag(A,-1)
C=diag(v,3)
B1 = tril(A,1)
B2 = triu(A,-1)
T1 =-triu(tril(ones(6,6),1),-1)+3*eye(6,6)
T2 =-diag(ones(5,1),-1) + diag(ones(6,1),0) - diag(ones(5,1),1)
T3 = toeplitz([2;-1;zeros(4,1)])
%-----------------------------------------------------------------------
Block Matrices
%-----------------------------------------------------------------------
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 
%-----------------------------------------------------------------------
  Review Cells
%-----------------------------------------------------------------------
% Script: cells
NYC_Lat  = struct('d',40,'m',45,'s',27);
NYC_Long = struct('d',75,'m',12,'s',32);
C1       = struct('name','New York','lat',NYC_Lat,'long',NYC_Long);
C1.name
C1.long.m
C=cell(2,2);
C{1,1}=[1 2;3 4];
C{1,2}=[5;6];
C{2,1}=[7 8];
C{2,2}=9;
M=[C{1,1} C{1,2};C{2,1} C{2,2}]
%-----------------------------------------------------------------------
% (1) matrix * vector  Sum A(:,i)*x(i)
Matrix*vector
%-----------------------------------------------------------------------
%     MatVecCo
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
%
% /$\       / --$----- \      -
% |$|       | --$----- |      -
% |$|   ___ | --$----- |      $
% |$|   \   | --$----- |      -
% |$| =  >  | --$----- |  *   -
% |$|   /__ | --$----- |      -
% |$|       | --$----- |      -
% |$|       | --$----- |      -
% \$/       \ --$----- / 
%-----------------------------------------------------------------------
%   matrix*vector  by rows A(i,:)*x  
%        MatVecRo   
%-----------------------------------------------------------------------
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
%       / -------- \     /$\
%       | -------- |     |$|
%   $ = | $$$$$$$$ |     |$|
%       | -------- |  *  |$|
%       | -------- |     |$|
%       | -------- |     |$|
%       \ -------- /     \$/ 
%-----------------------------------------------------------------------
