function varargout = irblsvds(varargin)
% IRBLSVDS: Finds a few singular values and singular vectors of a m x n matrix A.
%
% IRBLSVDS will find a few singular values and singular vectors such that A*V = U*S,
% where V is an orthonormal n x k matrix of "right" singular vectors, U is an
% orthonormal m x k matrix of "left" singular vectors, and S is a k x k diagonal 
% matrix of singular values. 
%
% IRBLSVDS uses the (m+n) x (m+n) Hermitian matrix  Z = [0  A; A' 0] and calls  
% IRBLEIGS(Z,OPTIONS) to find a few eigenvalues and eigenvectors of the matrix Z. 
% The singular values of A are the positive eigenvalues of Z, the "right" singular 
% vectors V, correspond to the last n elements of the eigenvectors of Z, and the 
% "left" singular vectors U, correspond to the first m elements of the eigenvectors of Z.
%
% [U,D,V,FLAG] = IRBLSVDS(A,OPTIONS) 
% or
% [U,D,V,FLAG] = IRBLSVDS('Afunc',M,N,OPTIONS) 
%
% The first input argument into IRBLSVDS can be a numeric matrix A or an M-file 
% ('Afunc') that  computes the product [0 A; A' 0] * X where X is a ((m+n) x blsz)
% matrix. If A is passed as a M-file then the first input argument is X, the second 
% input argument is M, the third input argument is N (A is a M x N matrix), and the 
% fourth input argument is blsz, i.e. Y = Afunc(X,m,n,blsz) or Y=Afun(X,m,n,blsz,FUNPAR). 
% See IRBLEIGS for a description of FUNPAR. In all the implementations IRBLSVDS(A, ...) 
% can be replaced with IRBLSVDS('Afunc',M,N,...).
%
% OUTPUT OPTIONS:
% ---------------
%
% I.)   IRBLSVDS(A) 
%       Displays the desired singular values.    
%
% II.)  S = IRBLSVDS(A)
%       Returns the desired singular values in the vector D. 
%
% III.) [U,S,V] = IRBLSVDS(A)
%       S is a diagonal matrix that contains the desired singular values in descending
%       order along the diagonal, the matrix V contains the corresponding "right" singular
%       vectors, and U contains the corresponding "left" singular vectors such that 
%       that A*V = U*S, V'*V = I, and U'*U = I.  
%
% IV.)  [U,S,V,FLAG] = IRBLSVDS(A) 
%       This option returns the same as (III) plus a variable FLAG that reports if 
%       the algorithm converges. If FLAG = 0 this implies normal return, 
%       all singular values have converged. If FLAG = 1 then the maximum number of 
%       iterations have been reached before all desired singular values have converged. 
%       
% INPUT OPTIONS:
% --------------
%                                   
%       ... = IRBLSVDS(A,OPTS) 
%       OPTS is a structure containing input parameters.
%
%  INPUT PARAMETER      DESCRIPTION 
%
%  OPTS.K            Number of desired singular values.             
%                    DEFAULT VALUE  K = 3
%
%  OPTS.SIGMA        Two letter string or numeric value specifying the location 
%                    of the desired singular values.            
%                    'LS'  Largest singular values.                 
%                    NVAL  A numeric value (>= 0). This routine searches for 
%                          singular values closest to the numeric value NVAL.
%                          Setting NVAL = 0 will search for the smallest
%                          singular values.
%                    DEFAULT VALUE   SIGMA = 'LS'
%
%  OPTS.TOL          Tolerance used for convergence. Convergence is determined when             
%                    || A*V - U*D ||_2 <= TOL*MAX(S). MAX(S) is approximated by
%                    largest absolute Ritz value of Z=[0  A; A' 0].  
%                    DEFAULT VALUE    TOL = 1d-6 
%
% The OPTS structure may also contain any of the following input parameters
% listed for IRBLEIGS;
%                     BLSZ, DISPR, ENDPT, FUNPAR, MAXIT, MAXDPOL, NBLS, SIZINT, 
%                     V0, and ZERTYP.
%
% See IRBLEIGS for more information.    
%
%  DATE:  09/16/02
%  VER:  1.0

%  AUTHOR:
%  James Baglama     University of Rhode Island, E-mail: jbaglama@math.uri.edu

% IRBLEIGS is required for IRBLSVDS to run. Verify that irbleigs is in the MATLAB path.
if exist('irbleigs') ~= 2, error('ERROR: The function IRBLEIGS is required for IRBLSVDS to run.'); end

% Persistent is used to keep A, BLSZ, M, and N in memory while IRBLSVDS calls
% itself to compute the matrix-vector products. The values are removed by the
% clear function before IRBLSVDS exits.
persistent A blsz m n;

% Incorrect number of output arguments requested.
if (nargout >= 5 | nargout==2 ) error('ERROR: Incorrect number of output arguments.'); end

%----------------------------%
% BEGIN: PARSE INPUT VALUES. %
%----------------------------%

% No input arguments, return help.
if nargin==0, help irblsvds, return, end

% Check input arguments and determine if IRBLSVDS is calling itself
% to compute the matrix-vector product. This avoids creating the
% matrix [0 A; A' 0] and any unnecessary products. This also allows IRBLEIGS
% to call the user input function 'Afunc' with input arguments 
% (X,m,n,blsz).
if nargin == 3 | nargin == 4
   if isnumeric(varargin{1}) & isnumeric(varargin{2}) & isnumeric(varargin{3})
      if (size(varargin{1},1) == n+m & size(varargin{1},2) == blsz & ...
         varargin{2} == n+m & varargin{3} == blsz)
         if isnumeric(A)
            varargout{1} = [A*varargin{1}(m+1:m+n,:);(varargin{1}(1:m,:)'*A)'];
            return;
         else
            if nargin == 3
               varargout{1} = feval(A,varargin{1},m,n,blsz);
            else
               varargout{1} = feval(A,varargin{1},m,n,blsz,varargin{4});
            end   
            return;    
         end  
      end   
   end
end   

%----------------------------%
% BEGIN: PARSE INPUT VALUES. %
%----------------------------%

% Get matrix A. Check type (numeric or character) and dimensions.
if (isstruct(varargin{1})), error('A must be a matrix.'), end
A = varargin{1};
if ischar(A)
   if nargin == 1, error('Need (M) for matrix A).'), end  
   m = varargin{2};
   if ~isnumeric(m) | length(m) ~= 1 
      error('Second argument M must be a numeric value.'); 
   end
   if nargin == 2, error('Need (N) for matrix A).'); end  
   n = varargin{3}; 
   if ~isnumeric(n) | length(n) ~= 1 
      error('Third argument N must be a numeric value.'); 
   end
else
   [m,n] = size(A);
   if ~isnumeric(A)
      error('A must be a numeric matrix.');
   end   
   if nnz(A) == 0 
      error('Matrix A contains all zeros.');
   end    
end

% Set K, MAXIT, TOL, SIGMA and BLSZ. 
K = 3; maxit = 2000; tol = 1d-6; sigma = 'LS';  
IK=[]; IS=[]; IB=[]; IM=[]; IT=[]; blsz = 3; 

% Get input arguments and adjust input values accordingly.
if nargin > 1 + 2*ischar(A) 
   OPTS = varargin{2+2*ischar(A):nargin};
   names = fieldnames(OPTS);
   IK = strmatch('K',upper(names),'exact');
   if  ~isempty(IK), K = getfield(OPTS,names{IK}); end 
   IB = strmatch('BLSZ',upper(names),'exact');
   if ~isempty(IB), blsz = getfield(OPTS,names{IB}); end 
   IS = strmatch('SIGMA',upper(names),'exact');
   if  ~isempty(IS), sigma = upper(getfield(OPTS,names{IS})); end
   IM = strmatch('MAXIT',upper(names),'exact');
   IT = strmatch('TOL',upper(names),'exact');
   if ~isempty(IT), tol = upper(getfield(OPTS,names{IT})); end
   I = strmatch('CHOLM',upper(names),'exact');
   if ~isempty(I), error('This option is not supported in IRBLSVDS.'), end
   I = strmatch('MPERM',upper(names),'exact');
   if ~isempty(I), error('This option is not supported in IRBLSVDS.'), end
   I = strmatch('EIGVEC',upper(names),'exact');
   if ~isempty(I), error('This option is not supported in IRBLSVDS.'), end
end  

% Check for input errors. 
if (~ischar(sigma) & ~isnumeric(sigma)), error('Unknown value for SIGMA.'), end
if K <= 0,       error('K must be a positive value.');    end   
if blsz <= 0,    error('BLSZ must be a positive value.'); end
if tol < 0,      error('TOL must be non-negative.');      end
if maxit <=0,    error('MAXIT must be positive.');        end
if K > min(n,m), K = min(n,m); end
        
% Modify SIGMA to be an acceptable input into IRBLEIGS.
if ischar(sigma)
   if length(sigma) ~= 2, error('Unknown value for SIGMA.'); end
   if strcmp(sigma,'LS')
      if ~isempty(IS), OPTS = setfield(OPTS,names{IS},'LE'); else OPTS.sigma='LE'; end
   else
      error('Unknown value for SIGMA.')
   end   
else   
   % Set K = BLSZ+2*K to help ensure K singular values are found when SIGMA is numeric.
   if sigma < 0, OPTS = setfield(OPTS,names{IS},0); end
   if ~isempty(IK),OPTS = setfield(OPTS,names{IK},2*K+blsz); else OPTS.K= 2*K+blsz; end
end

% Set default value for MAXIT to be 300, if MAXIT is not part of the structure.
if isempty(IM), OPTS.maxit=maxit; end

%--------------------------%
% END: PARSE INPUT VALUES. %
%--------------------------%

%----------------------------------------------------%
% BEGIN: CALL IRBLEIGS, SORT RESULTS, OUTPUT VALUES. %
%----------------------------------------------------%

% Rmax is used in IRBLEIGS and holds the maximum absolute value of 
% all computed Ritz values of the matrix Z = [0  A; A' 0] and is used 
% to approximate the largest singular value of A.           
global Rmax; Rmax=[];

% Call IRBLEIGS to compute the eigenvalues of the matrix [0 A; A' 0].
[X,S,FLAG] = irbleigs('irblsvds',n+m,OPTS); 

% Extract the singular values from S. Sort the values, remove all negative
% eigenvalues, and determine which zero eigenvalues are zero singular values.
[U,S,V] = extsingvals(K,m,n,Rmax,sigma,S,X,tol);
if length(S) < K | FLAG(1)~=0, 
    disp('WARNING: One or all of the following may have occurred.');
    disp('     1.) Not all singular values converged.');
    disp('     2.) Not all desired singular values were found.');
    disp('     3.) Singular value incorrectly labeled largest, smallest or closest to SIGMA.');
    FLAG(1) = 1;
end 
  
% Output results.
if nargout == 0, SingularValues = S, end
if nargout == 1, varargout{1} = S;   end
if nargout > 1
   varargout{1} = U; varargout{2} = diag(S); varargout{3}=V;
   if nargout == 4, varargout{4}=FLAG(1); end
end 

% Reinitializes persistent variables by clearing the function IRBLSVDS.
clear irblsvds;

%--------------------------------------------------%
% END: CALL IRBLEIGS, SORT RESULTS, OUTPUT VALUES. %
%--------------------------------------------------%
 
%--------------------%
% BEGIN: EXTSINGVALS %
%--------------------%

function [U,S,V] = extsingvals(K,m,n,Rmax,sigma,S,X,tol)
% This function extracts the singular values from S. Sort the values, removes all negative
% eigenvalues, determines if any zero eigenvalues are singular values, and outputs the
% singular vectors.

% James Baglama
% 09/16/02

% Initialization and description of local variables.
JPOS = [];  % Variable that holds which eigenvalues of Z are considered positive and hence singular values of A.
JZER = [];  % Variable that holds which eigenvalues of Z are considered zero.

% Extract the diagonal elements from S;
S=diag(S);

% Determine the positive singular values.
JPOS = find(S > sqrt(eps)*max(n,m)*Rmax);

% Determine the zero eigenvalues of Z.  
JZER = find(abs(S) <= sqrt(eps)*max(n,m)*Rmax);

% Extract vectors U and V from X. The vectors U and V have 2-norm equal to 1/sqrt(2).
V = sqrt(2)*X(m+1:n+m,:); U = sqrt(2)*X(1:m,:);
   
% Determine which zero eigenvalues of [0 A; A' 0] are also zero 
% singular values. By construction the matrix [0 A; A' 0] will 
% have abs(m-n) zero eigenvalues that are not zero singular values.
%
% !!! WARNING: This test is not fail safe !!!!
%
% A spurious zero singular value can result from this test.
if ~isempty(JZER)
   if m >= n
      if rank(V(:,JZER),sqrt(eps)*Rmax) > 0
         [Q,R] = qr(V(:,JZER),0); 
         JZER = JZER(find(abs(diag(R)) >= max(sqrt(eps),tol)*max(n,m)*Rmax));
      else
         JZER=[];  
      end
   else
      if rank(U(:,JZER),sqrt(eps)*Rmax) > 0
         [Q,R] = qr(U(:,JZER),0);
         JZER = JZER(find(abs(diag(R)) >= max(sqrt(eps),tol)*max(n,m)*Rmax));
      else
         JZER=[];  
      end 
   end 
end  
V(:,JZER) = orth(V(:,JZER)); U(:,JZER) = orth(U(:,JZER));
J = union(JPOS,JZER); S = abs(S(J)); V = V(:,J); U = U(:,J);
if ischar(sigma)
    [sortval,I] = sort(S); I = flipud(I); I=I(1:min(K,length(S)));
else
    [sortval,I] = sort(abs(S-sigma)); I=I(1:min(K,length(S)));
end 
[sortval,J] = sort(S(I)); I = I(J); S = S(I); V = V(:,I); U = U(:,I);

%------------------%
% END: EXTSINGVALS %
%------------------%