function [ Vectors, Values, Psi ] = pc_vectors ( A, numvecs )

%*****************************************************************************80
%
%% PC_VECTORS gets principal component vectors.
%
%  Discussion:
%
%    Get the top numvecs eigenvectors of the covariance matrix
%    of A, using Turk and Pentland's trick for numrows >> numcols
%    Returns the eigenvectors as the colums of Vectors and a
%    vector of ALL the eigenvalues in Values.
%
%  Modified:
%
%    03 December 2018
%
%  Modifier:
%
%    John Burkardt
%
%  Author:
%
%    Matthew Dailey 2000
%
%  Parameters:
%
%    Input, real A(M,N), an array of data.
%
%    Input, integer NUMVECS, the number of vectors requested.
%    NUMVECS <= N.
%
%    Output, real VECTORS(M,NUMVECS), the requested vectors.
%
%    Output, real VALUES(N), all the eigenvalues of the covariance matrix.
%
%    Output, real PSI(M), the average vector.
%
  [ m, n ] = size ( A );
%
%  Compute the average.
%
  Psi = mean ( A' )';
%
%  Subtract the average from each column of A.
%
  for j = 1 : n
    A(:,j) = A(:,j) - Psi;
  end
%
%  Compute the covariance matrix.
%
  L = A' * A;
%
%  Eigenvectors and eigenvalues of covariance matrix.
%
  [ Vectors, Values ] = eig ( L );
%
%  Sort the values according to size, and sort the vectors accordingly.
%
  Values = diag ( Values );
  Values = Values / ( n - 1 );
  [ Values, sort_index ] = sort ( Values, 'descend' );

  Vectors = Vectors(:,sort_index);

% [Vectors,Values] = sortem2(Vectors,Values);
%
%  Convert the eigenvectors of A'*A into eigenvectors of A*A'
%
  Vectors = A * Vectors;
%
%  Normalize Vectors.
%
  for j = 1 : n
    Vectors(:,j) = Vectors(:,j) / norm ( Vectors(:,j) );
  end
%
%  Kill vectors corr. to tiny evalues
%
  num_good = 0;

  for j = 1 : n
    if ( abs ( Values(j) ) < 0.00001 )
      Values(j) = 0.0;
      Vectors(:,j) = zeros ( m, 1 );
    else
      num_good = num_good + 1;
    end
  end

  if ( num_good < numvecs )
    fprintf(1,'Warning: numvecs is %d; only %d exist.\n',numvecs,num_good);
    numvecs = num_good;
  end;

  Vectors = Vectors(:,1:numvecs);

  return
end
