function r = hankel_cholesky_upper ( n, h )

%*****************************************************************************80
%
%% hankel_cholesky_upper() computes the upper Cholesky factor of a Hankel matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 January 2017
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Phillips,
%    The triangular decomposition of Hankel matrices,
%    Mathematics of Computation,
%    Volume 25, Number 115, July 1971, pages 599-602.
%
%  Input:
%
%    integer N, the order of the matrix.
%
%    real H(2*N-1), the values defining the antidiagonals.
%
%  Output:
%
%    real R(N,N), the upper triangular Cholesky factor.
%
  h = h(:);
  a = zeros ( n, 1 );
  b = zeros ( n, 1 );
  c = zeros ( 2 * n - 1, 2 * n - 1 );

  c(1,1:2*n-1) = h(1:2*n-1)';

  for i = 1 : n - 1

    if ( i == 1 )
      a = c(1,2) / c(1,1);
      b = 0.0;
    else
      a = c(i,i+1) / c(i,i) - c(i-1,i) / c(i-1,i-1);
      b = c(i,i) / c(i-1,i-1);
    end

    for j = i + 1 : 2 * n - i - 1
      c(i+1,j) = c(i,j+1) - a * c(i,j);
      if ( 1 < i )
        c(i+1,j) = c(i+1,j) - b * c(i-1,j);
      end
    end
  end

  r(1:n,1:n) = c(1:n,1:n);
%
%  Normalize.
%
  for i = 1 : n
    t = sqrt ( r(i,i) );
    r(i,1:i-1) = 0.0;
    r(i,i:n) = r(i,i:n) / t;
  end

  return
end

