function Glambda=gcvgen(lambda,n,s2,t,sumt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Glambda=gcvgen(lambda,n,s2,t,sumt);
% evaluates crossvalidation function at value lambda for scalar regularization parameter 
%
%       G(lambda)= || P(lambda) Y ||^2 ./ Trace(P(lambda))^2 
%
% for use in Tikhonov regularization
%
%         min {|| A x - Y ||^2+lambda^2 || x ||^2} (function tikh), 
%
% where        
%
%      P(lambda) = I - A inv(A' A + lambda^2 I) A'
%
% and
%
%     A = U diag(s) V' by compact SVD (csvd). 
%
% The seudoinverse of A' A, A' A+lambda^2 I, is implicitly used.
% 
% Input arguments:
%
%   lambda  - optimal regularization parameter
%   n       - number of rows of Y
%   s2      - column vector of squares of singular values of A truncated for numerical stability
%   t       - sum of squares of rows of the product of transposed left orthogonal matrix from 
%           - SVD of A by data matrix Y
%
% Output argument:
%
%   Glambda - G(lambda)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
lambda = 10.^lambda;
fi=s2./(s2+lambda.^2);
num=sumt+t'*((1-fi).^2);
den=(n-sum(fi)).^2;
if den < num*eps
  Glambda = realmax;
else
  Glambda=num/den;
end
