function splus = compute_splus ( s )

%% COMPUTE_SPLUS computes S+, the "inverse" of the S factor.
%
%  Discussion:
%
%    S is an MxN diagonal matrix, whose diagonal entries are nonnegative,
%    and which appear on the diagonal in descending order, with zero values
%    last.
%
%    SPLUS is an NxM diagonal matrix.  Each nonzero entry SPLUS(I,I) is 
%    the inverse of the corresponding nonzero entry S(I,I).
%
%    If A has the SVD factorization A = U * S * V', then the pseudoinverse
%    APLUS is V * SPLUS * U'.
%
%    Note that the nonzero diagonal values in SPLUS appear in ASCENDING order,
%    so that, strictly speaking, V * SPLUS * U' is not a singular value
%    decomposition. 
%
%    We include a tolerance TOL which allows us to ignore relatively
%    small entries of S.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(M,N), the singular value matrix.
%
%    Output, real SPLUS(N,M), the pseudoinverse of S.
%

%
%  Get the size of S.
%
  [ m, n ] = size ( s );
%
%  Think of SPLUS as a row vector.
%
  splus = reshape ( s, 1, m * n );
%
%  Locate the nonzero elements in the row vector.
%
  tol = s(1,1) * sqrt ( eps );

  nonzeros = find ( splus > tol );
%
%  Replace the nonzero elements by their inverse.
% 
  splus(nonzeros) = 1.0 ./ splus(nonzeros);
%
%  Think of the row vector as an MxN matrix.
%
  splus = reshape ( splus, m, n );
%
%  Tranpose the matrix.
%
  splus = splus';

  return
end

  
