function [s,v,u] = OperatorNorm(x,tol)

% Input:    x       initial right singular vector
%           tol     desired absolute tolerance
% Output:   s       maximal singular value (operator norm)
%           v       right singular vector
%           u       left  singular vector

n = length(x); y = zeros(n,1);
k = 1:n;
lambda = 0; lambda0 = 1; 
while abs(lambda-lambda0) > tol
    xhat = x/norm(x);
    for j=1:n, y(j)=(1./((j+k-1).*(j+k)/2-(k-1)))*xhat; end
    for j=1:n, x(j)=(1./((k+j-1).*(k+j)/2-(j-1)))*y;    end
    lambda0 = lambda; lambda   = x'*xhat;
end
s = sqrt(lambda);
v = x/norm(x);
u = y/norm(y);

return