function [V, d, r] = peigs(A, rmax)%PEIGS   Finds positive eigenvalues and corresponding eigenvectors. %  %    [V, d, r] = PEIGS(A, rmax) determines the positive eigenvalues d%    and corresponding eigenvectors V of a matrix A. The input%    parameter rmax is an upper bound on the number of positive%    eigenvalues of A, and the output r is an estimate of the actual%    number of positive eigenvalues of A. The eigenvalues are returned%    as the vector d(1:r). The eigenvectors are returned as the%    columns of the matrix V(:, 1:r).%%    PEIGS calls the function EIGS to compute the first rmax%    eigenpairs of A by Arnoldi iterations. Eigenpairs corresponding%    eigenvalues that are nearly zero or less than zero are%    subsequently discarded.%%    PEIGS is only efficient if A is strongly rank-deficient with%    only a small number of positive eigenvalues, so that rmax <<%    min(size(A)). If A has full rank, EIG should be used in place of%    PEIGS.%%    See also: EIGS.  error(nargchk(2, 2, nargin))                    % check number of input arguments   [m, n]  = size(A);   if rmax > min(m,n)    rmax  = min(m,n);                             % rank cannot exceed size of A   end    % get first rmax eigenvectors of A  warning off  eigsopt.disp = 0;                               % do not display eigenvalues  [V, d]       = eigs(A, rmax, 'lm', eigsopt);  warning on    % output of eigs differs in different Matlab versions  if prod(size(d)) > rmax    d          = diag(d);                         % ensure d is vector  end  % ensure that eigenvalues are monotonically decreasing  [d, I]       = sort(d, 'descend');  V            = V(:, I);    % estimate number of positive eigenvalues of A  d_min        = max(d) * max(m,n) * eps;   r            = sum(d > d_min);				   % discard eigenpairs with eigenvalues that are close to or less than zero  d            = d(1:r);  V            = V(:, 1:r);  d            = d(:);				  % ensure d is column vector