function x = BlkTriSolve(A,C,D,b)
  % A is length-p cell array with entries that are r-by-r matrices.
  % C and D are length-(p-1) cell arrays with entries that are r-by-r matrices
  % b is an (pr)-by-1 vector.
  % x solves Tx = b where T is block tridiagonal with diagonal blocks 
  %     A{1},...,A{p}, subdiagonal blocks C{1},...C{p-1}, and superdiagonal blocks
  %     D{1},...,D{p-1}.
  p = length(A);
  [r,r] = size(A{1});
  n = p*r;
  E = cell(p-2,1);
  for k=1:p-1
      % Compute QR factorization of the 2r-by-r matrix
      M1 = [A{k};C{k}];
      [Qk,Rk] = qr(M1);
      
      % Update A{k} and b...
      A{k} = Rk(1:r,1:r);
      rows = 1+(k-1)*r:(k+1)*r;
      b(rows) = Qk'*b(rows);
      
      % Perform the update Qk'*M2 where
      M2 = [D{k};A{k+1}];
      M2 = Qk'*M2; 
      D{k}   = M2(1:r,:); 
      A{k+1} = M2(r+1:2*r,:);
      
      if k<=p-2
         % Exploiting zeros, perform the update 
         %       E{k}             zeros(r,r)
         %               =  Qk' *              
         %       D{k+1}           D{k+1}
         M3 = Qk(r+1:2*r,:)'*D{k+1};
         E{k} = M3(1:r,:); D{k+1} = M3(r+1:2*r,:);
      end
  end
  % Process the last diagonal block...
  [Qk,Rk] = qr(A{p});
  A{p} = Rk; 
  b(n-r+1:n) = Qk'* b(n-r+1:n);
  
  % Now back solve
  x = zeros(n,1);
  for k=p:-1:1
      rows = 1+(k-1)*r:k*r;
      if k<=p-1
          % This block row has a D-block...
          b(rows) = b(rows) - D{k}*x(r+rows);
      end
      if k<=p-2
          % This block row has an E-block
          b(rows) = b(rows) - E{k}*x(2*r+rows);
      end
      % A{k} is upper triangular so...
      x(rows) = A{k}\b(rows);
  end
  
  % -1 if do not fully exploit zeros
  % -3 if you form T as an array