function [gamma, delta, p, q] = lbidiag2 (G, d, gamma, delta, p, q, k)

% ------------------------------------------------------------------------
% Lanczos Bidiagonalization II Algorithm
%
% References:
% [1] G. H. GOLUB and U. von MATT,
%     "Generalized Cross-Validation for Large-Scale Problems", 1996,
%     TR-96-28
%
% [2] G. H. GOLUB and U. von MATT,
%     "Tikhonov Regularization for Large Scale Problems", 1997, SCCM-97-03
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Initialization
% ------------------------------------------------------------------------
[m,n] = size(G);                % dimension of matrix G
normG = normest(G);             % estimated norm of G
kmax  = sqrt(min(m,n));         % upperbound on number of iterations
kmin  = ceil(3*log(min(m,n)));  % lowerbound on number of iterations
tol   = max(m,n)*10*eps*normG;  % stopping tolerance

if nargin <= 2
  gamma = [];                   % result: above diagonal
  delta = [];                   % result: diagonal
  k     = 1;                    % iteration index
  p     = d/norm(d, 2);         % starting vector
else
  kmax  = 2*k;
end

% ------------------------------------------------------------------------
% Loop
% ------------------------------------------------------------------------
while k <= kmax % TODO:!!!

  if (k <= 1), q = G'*p;
  else,        q = G'*p - delta(k-1)*q; end

  gamma(k) = norm(q, 2);
  
  % termination criteria
  if k > kmin && abs (gamma(k)) <= tol, gamma = gamma(1:k-1); break; end

  q        = q / gamma(k);
  p        = G*q - gamma(k)*p;
  delta(k) = norm(p, 2);

  % termination criteria
  if k > kmin && abs (delta(k)) <= tol, break; end

  p        = p / delta(k);
  k        = k+1;

end
