  function P1Grade()
% Checks out HermHactor
% Generate a random nxn tridiagonal hermitian positive definite matrix A
% and represent it with three real vectors
n = 100000;

% Generate D and L
d = rand(n,1);
c = randn(n-1,1);
e = randn(n-1,1);

%              1           0            0          0
%  L  =    c(1)+ie(1)      1            0          0
%              0       c(2)+ie(2)       1          0
%              0           0        c(3)+ie(3)     1
%
%             d(1)         0            0          0
%  D  =        0          d(2)          0          0
%              0           0           d(3)        0
%              0           0            0         d(4)
%
%              1       c(1)-ie(1)            0          0
%  L'  =       0           1            c(2)-ie(2)      0
%              0           0                 1      c(3)-ie(3)
%              0           0                 0          1
%  
%           a(1)     g(1)-ih(1)         0              0
%   A =  g(1)+ih(1)     a(2)        g(2)-ih(2)         0
%            0       g(2)+ih(2)        a(3)        g(3)-ih(3)
%            0           0          g(3)+ih(3)         a(4)
%
% Comparing coefficients in A = LDL'
%
%        a(1)        = d(1)
%        g(1)+ih(1)  = d(1)c(1)        + id(1)e(1)
%        a(2)        = d(2)+d(1)c(1)^2 + d(1)e(1)^2
%        g(2)+ih(2)  = d(2)c(2)        + id(2)e(2)
%        a(3)        = d(3)+d(2)c(2)^2 + d(2)e(2)^2
%        g(3)+ih(3)  = d(3)c(3)        + id(3)e(3)
%        a(4)        = d(4)+d(3)c(3)^2 + d(3)e(3)^2
%
% Generate A, i.e.,
  a = zeros(n,1);
  g = zeros(n-1,1);
  h = zeros(n-1,1);
  a(1) = d(1);
  for k=1:n-1
      g(k) = d(k)*c(k);
      h(k) = d(k)*e(k);
      a(k+1) = d(k+1) + d(k)*(c(k)^2+e(k)^2);
  end
  % Generate L and D is sparse format...
  [d0,c0,e0] = HermFactor(a,g,h);
  L = spdiags([[c0;0]+sqrt(-1)*[e0;0] ones(n,1)],[-1 0],n,n);
  D = spdiags(d0,0,n,n);
  % Generate A in sparse format
  A = spdiags([[g;0]+sqrt(-1)*[h;0] a [0;g]-sqrt(-1)*[0;h]],[-1 0 1],n,n);
  clc 
  % Error...
  fprintf('|| A - LDL'' || = %10.2e\n',norm(A-L*D*L',1))
  % Benchmark...
  nRepeat = 200;
  tYOU = 0;
  tic
  for k=1:nRepeat
     [d0,c0,e0] = HermFactor(a,g,h);
  end
  tYOU = toc;
  
  tCVL = 0;
  tic
  for k=1:nRepeat
     [d0,c0,e0] = HermFactorCVL(a,g,h);
  end
  tCVL = toc;
  fprintf('\nn = %1d\n',n)
  fprintf('tYOU = %5.3f\n',tYOU)
  fprintf('tCVL = %5.3f\n',tCVL)
  




  function [d,c,e] = HermFactorCVL(a,g,h)
% a is a real column n-vector.
% g and h are real column (n-1)-vectors.
% Let A be the n-by-n Hermitian tridiagonal matrix with diagonal a and subdiagonal 
% g + i*h,  where i^2 = -1. Assume that A is positive definite so that it has
% the factorization A = L*D*L' where D is diagonal and L is unit lower bidiagonal.
% 
% d is a real column n-vector so D = diag(d).
% c and e are real column (n-1)-vectors so that the subdiagonal
% of L is  c + i*e.

% Using the above n=4 recipes, here is the algorithm...
%        d(1) = a(1)
%        c(1) = g(1)/d(1)   
%        e(1) = h(1)/d(1)
%        d(2) = a(2) -d(1)c(1)^2 - d(1)e(1)^2 = a(2) - g(1)c(1) - h(1)e(1)
%        c(2) = g(2)/d(2)   
%        e(2) = h(2)/d(2)
%        d(3) = a(3) -d(2)c(2)^2 - d(2)e(2)^2 = a(3) - g(2)c(2) - h(2)e(2)
%        c(3) = g(3)/d(3)   
%        e(3) = h(3)/d(3)
%        d(4) = a(4) -d(3)c(3)^2 - d(3)e(3)^2 = a(4) - g(3)c(3) - h(3)e(3)

n = length(a);
d = zeros(n,1);
c = zeros(n-1,1);
e = zeros(n-1,1);
d(1) = a(1);
for k=1:n-1
    c(k) = g(k)/d(k);  
    e(k) = h(k)/d(k);
    d(k+1) = a(k+1) - g(k)*c(k) - h(k)*e(k);
end


% -1, -2 for gross inefficiencies