  function x = BlockTriSolve(A,f,b)
% A is a p-by-1 cell array with A{i) a q-by-q symmetric matrix
% f is a (p-1)-by-1 array of real numbers
% Let M be a p-by-p block tridiagonal matrix with diagonal blocks
%                   M_{ii} = A{i}   i=1:p
% and  subdiagonal blocks
%                   M_{i+1,i} = f(i)ei*eq' i=1:p-1
% where e1 and eq are the first and last columns of eye(q,q).
% 
% The matrix M is positive definite.
% b is a p-by-1 cell array with b{i} a q-by-1 vector
% x is a p-by-1 cell array with x{i} a q-by-1 vector
% M[x{1};...;x{p}] = [b{1};...;b{p}]

  p = length(b);
  q = length(b{1});
  
% Let G be the cholesky factorization of M. G is block lower
% bidiagonal. Denote the diagonal blocks by G_{1},...,G_{p}
% and the subdiagonal blocks by H_{1},...,H_{p-1}. Comparing blocks we get 
%
%        G_{1}G_{1}' = A{1}  (Cholesky)
%        for i=1:p-1
%            Solve H_{i}G_{i}' = f(i)ei*eq' for H_{i}
%            G_{i+1}G_{i+1}' = A_{i+1} - H_{i}H_{i}'  (Cholesky)
%        end
%      
% It follows that H_{i} = (f(i)/G_{i}(q,q)) ei*eq' and so we obtain

  A{1} = chol(A{1})';
  for i=1:p-1
      f(i) = f(i)/A{i}(q,q);
      A{i+1}(1,1) = A{i+1}(1,1) - f(i)^2;
      A{i+1} = chol(A{i+1})';
  end
  
% Block Forward elimination Gy = b..
  b{1} = A{1}\b{1};
  for i=1:p-1
      b{i+1}(1) = b{i+1}(1) - f(i)*b{i}(q);
      b{i+1} = A{i+1}\b{i+1};
  end
  
% Block Back Substitution G'x= y
  b{p} = A{p}'\b{p};
  for i=p-1:-1:1
      b{i}(q) = b{i}(q) - f(i)*b{i+1}(1);
      b{i} = A{i}'\b{i};
  end
  
  x = b;
  
  % No points deducted i you made f a cell array.
  % my version of matlab does not support chol(  , 'lower'). Ignored this.

