function [ nf, xf, wf, vf ] = line_fekete_monomial ( m, a, b, n, x )

%*****************************************************************************80
%
%% line_fekete_monomial() computes approximate Fekete points in an interval [A,B].
%
%  Discussion:
%
%    We use the uniform weight and the monomial basis:
%
%      P(j) = x^(j-1)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 October 2022
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alvise Sommariva, Marco Vianello,
%    Computing approximate Fekete points by QR factorizations of Vandermonde 
%    matrices,
%    Computers and Mathematics with Applications,
%    Volume 57, 2009, pages 1324-1336.
%    
%  Input:
%
%    integer M, the number of basis polynomials.
%
%    real A, B, the endpoints of the interval.
%
%    integer N, the number of sample points.
%    M <= N.
%
%    real X(N), the coordinates of the sample points.
%
%  Output:
%
%    integer NF, the number of Fekete points.
%    If the computation is successful, NF = M.
%
%    real XF(NF), the coordinates of the Fekete points.
%
%    real WF(NF), the weights of the Fekete points.
%
%    real VF(NF,M), the nonsingular Vandermonde submatrix.
%
  if ( n < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'line_fekete_monomial(): Fatal error!\n' );
    fprintf ( 1, '  N < M.\n' );
    error ( 'line_fekete_monomial(): Fatal error!' );
  end
%
%  Destroy all row vectors!
%
  x = x(:);
%
%  Moments of monomials:
%
  mom = line_monomial_moments ( a, b, m );
%
%  Form the rectangular Vandermonde matrix V for the polynomial basis.
%
  v = zeros ( m, n );

  v(1,1:n) = 1.0;
  for i = 2 : m
    v(i,1:n) = v(i-1,1:n) .* x(1:n)';
  end
%
%  Solve the system for the weights W.
%  Note that OCTAVE backslash does not work the way MATLAB's does.
%  So we have to use QR_FACTOR_SOLVE instead.
%
% w = v \ mom;

  w = qr_factor_solve ( m, n, v, mom );
%
%  Locate nonzero W's.
%
  ind = find ( w ~= 0.0 );
%
%  Retain data associated with nonzero W's.
%
  nf = length ( ind );
  xf = x ( ind );
  wf = w ( ind );
  vf = v ( :, ind );

  return
end

