function [U, S, V] = lyngby_svd(X, arg1, arg2, arg3, arg4, arg5, ...
    arg6, arg7, arg8)

% lyngby_svd           - Singular value decomposition
%
%       function [U, S, V] = lyngby_svd(X, ...
%          'PropertyName',PropertyValue)
%
%       Input:    X   Matrix. 
%
%       Property: ReduceComponents [ {all} | Integer ] Number of
%                                  components to reduce to
%
%       Output:   U   The eigensequence
%                 S   Singular values in a vector
%                 V   The eigenimages
%
%       Perform singular value decomposition: 
%
%                      X = U * diag(S) * V'
%
%       The returned number of components is the minimum of the rank
%       of X and the number specified with the property
%       'ReduceComponents'. 
%
%       The singular value decomposition is performed as an eigenvalue
%       decomposition on either the inner or the outer product (X*X'
%       or X'*X) depending on which one is the computational most
%       favorable. 
%
%       See also LYNGBY, LYNGBY_KM_MAIN, LYNGBY_SOP_MAIN,
%                LYNGBY_UI_SVD_INIT, LYNGBY_UIS_SVD. 
%
% $Id: lyngby_svd.m,v 1.11 2002/11/29 18:00:36 fnielsen Exp $


    % Default values
    reduceComponents = -1; 

    % Parse properties
    n = 1;
    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 
      if strcmp(arg, 'reducecomponents')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	reduceComponents = arg;
	if reduceComponents < 1
	  error('The number of ''ReduceComponents'' should be an integer larger than 1.')
	end
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    [rX, cX] = size(X);

    if rX <= cX

      lyngby_log('Beginning computation of eigenvalues ...')
      [U,S] = eig(X*X');
      lyngby_log('Computation of eigenvalues completed!');

      [S,I] = sort(-diag(S));
      S = -S;
      
      % Strip all null-definite directions
      tol = length(S) * S(1) * eps;
      posdef = max(find(S>tol));
      if reduceComponents ~= -1
	reduceComponents = min([reduceComponents posdef]);
      elseif posdef < length(S)
	reduceComponents = posdef;
      end
      
      if reduceComponents ~= -1
	reduceIndices = 1:reduceComponents;
	U = U(:,I(reduceIndices));
	S = S(reduceIndices, 1);
      else
	U = U(:,I); 
      end

      S = sqrt(S);
      V = X' * U * diag(1./S);
      
    else 

      lyngby_log('Beginning computation of eigenvalues ...')
      [V,S] = eig(X'*X);
      lyngby_log('Computation of eigenvalues completed!');

      [S,I] = sort(-diag(S));
      S = -S;
      
      % Strip all null-definite directions
      tol = length(S) * S(1) * eps;
      posdef = max(find(S>tol));
      if reduceComponents ~= -1
	reduceComponents = min([reduceComponents posdef]);
      elseif posdef < length(S)
	reduceComponents = posdef;
      end

      if reduceComponents ~= -1
	reduceIndices = 1:reduceComponents;
	V = V(:,I(reduceIndices));
	S = S(reduceIndices, 1);
      else
	V = V(:,I); 
      end

      S = sqrt(S);
      U = X * V * diag(1./S);

    end












