
    function [L,K] = SymSkew(A,C)
  % Computes the symskew factorization of a symskew matrix T encoded by A and C.
  % A is a length-p cell array where A{i} is an r-by-r symmetric positive
  %    definite matrix,i=1:p. These are the diagonal blocks.
  % C is a length-(p-1) cell array where C{i} is and r-by-r matrix, i=1:p-1.
  %    These are the subdiagonal blocks.
  % L is a length-p cell array where L{i} is r-by-r and lower triangular, i=1:p.
  % K is a length-(p-1) cell array where K{i} is r-by-r, i=1:p-1.
  % L and K encode the symskew factorization of T.
  
  p = length(A); L = cell(p,1); K = cell(p-1,1);
  % Comparing (1,1) blocks L{1}*L{1}' = A{1} so
  % L{1} = chol(A{1},'lower');
  L{1} = chol(A{1})';
  for i=1:p-1
      % Comparing (i+1,i) blocks says K{i}L{i}' = C{i}, 
      % i.e., L{i}K{i}' = C{i}'
      K{i} = (L{i}\C{i}')';
      % Comparing (i+1,i+1) blocks L{i+1}L{i+1}' - K{i}K{i}' = A{i+1}  so
      AiTilde = A{i+1}+K{i}*K{i}';
      %L{i+1} = chol(AiTilde,'lower');
      L{i+1} = chol(AiTilde)';
  end
  
  % -1 for deficient comments...
  
  % The reason why AiTilde is positive definite:
  
  % For any nonzero r-vector v we have
  %    v'*AiTilde*v = v'*(A{i+1} + K{i}*K{i}')*v 
  %                 = v'*A{i}*v  + v'*K{i}*K{i}*v
  %                 = v'*A{i}*v  + norm(K{i}'*v,2)^2
  %                 >= v'*A{i}*v 
  %                 > 0
  % The cannot have the Ai's too badly conditioned compared tothe Ci's