% % CS 4220 (Spring 2010) Assignment 3
% Test Script P1. Checks out BlockTriDiag and BlockSol.

clc
randn('seed',0)
rand('seed',0)

% Check BlockTriDiag on a small-N example...
N = 3; m = 10; A = cell(N,N); C = cell(N-1,N-1);
for i=1:N
    M = randn(m,m); A{i} = M'*M + 2*m*eye(m); 
    if i<N
      C{i} = rand(m,m);
    end
end
[G,H] = BlockTriDiag(A,C);
Z = zeros(m,m); 
% Does Atilde = Gtilde*Gtilde' where
Atilde = [ A{1} C{1}' Z ; C{1} A{2} C{2}' ; Z C{2} A{3}];
Gtilde = [ G{1}  Z    Z ; H{1} G{2}  Z    ; Z H{2} G{3}];
errBlockTriDiag = norm(Atilde - Gtilde*Gtilde')

% Check BlockTriDiagSol...
N = 100; m = 50; A = cell(N,1); C = cell(N-1,1); b = cell(N,1);
for i=1:N
    M = randn(m,m); A{i} = M'*M + 2*m*eye(m);
    b{i} = randn(m,1);
    if i<N
      C{i} = rand(m,m);
    end
end
[G,H] = BlockTriDiag(A,C);
x = BlockSol(G,H,b);
% How small is Atilde*x-b?
res = norm(A{1}*x{1} + C{1}'*x{2} - b{1});
for i=2:N-1
    res = res + norm(C{i-1}*x{i-1} + A{i}*x{i} + C{i}'*x{i+1} - b{i});
end
res = res + norm(C{N-1}*x{N-1} + A{N}*x{N} - b{N})


