function l = hankel_spd_cholesky_lower ( n, lii, liim1 )

%*****************************************************************************80
%
%% hankel_spd_cholesky_lower() returns L such that L*L' is Hankel SPD.
%
%  Discussion:
%
%    In other words, H = L * L' is a symmetric positive definite (SPD) matrix
%    with the property that H is constant along antidiagonals, so that
%
%      H(I+J) = h(k-1), for 1 <= I, J <= N, 1 <= K <= 2*N-1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 January 2017
%
%  Author:
%
%    S Al-Homidan, M Alshahrani.
%    MATLAB implementation by John Burkardt.
%
%  Reference:
%
%    S Al-Homidan, M Alshahrani,
%    Positive Definite Hankel Matrices Using Cholesky Factorization,
%    Computational Methods in Applied Mathematics,
%    Volume 9, Number 3, pages 221-225, 2009.
%
%  Input:
%
%    integer N, the order of the matrix.
%
%    real LII(N), values to be used in L(I,I), for 1 <= I <= N.
%
%    real LIIM1(N-1), values to be used in L(I+1,I) for 1 <= I <= N-1.
%
%  Output:
%
%    real L(N,N), the lower Cholesky factor.
%
  l = zeros ( n, n );

  for i = 1 : n
    l(i,i) = lii(i);
  end

  for i = 1 : n - 1
    l(i+1,i) = liim1(i);
  end

  for i = 3 : n
    for j = 1 : i - 2

      if ( mod ( i + j, 2 ) == 0 )
        q = ( i + j ) / 2;
        r = q;
      else
        q = ( i + j - 1 ) / 2;
        r = q + 1;
      end

      alpha = 0.0;
      for s = 1 : q
        alpha = alpha + l(q,s) * l(r,s);
      end

      beta = 0.0;
      for t = 1 : j - 1
        beta = beta + l(i,t) * l(j,t);
      end

      l(i,j) = ( alpha - beta ) / l(j,j);

    end
  end

  return
end

