function varargout = irlbablk(varargin)
%
% IRLBABLK will find a few singular values and singular vectors of A 
% such that A*V = U*S and A'*U = V*S, where V is an orthonormal n x k matrix 
% of "right" singular vectors, U is an orthonormal m x k matrix of "left" singular 
% vectors, and S is a k x k diagonal matrix of singular values. This program 
% utilizes the a Augmented Block Lanczos bidiagonalization method to find the
% singular triples. The program uses only matrix vector products to find the
% singular triples at either end of the spectrum. 
%
% [U,S,V,FLAG] = IRLBABLK(A,OPTIONS) or
% [U,S,V,FLAG] = IRLBABLK('Afunc',M,N,OPTIONS) 
%
% The first input argument into IRLBABLK can be a numeric matrix A or an M-file 
% ('Afunc'). If the M x N matrix A is a filename then Y = Afunc(X,M,N,BLSZ,'transpose'). 
% If transpose = 'F', then X is an (N x BLSZ) matrix and Y = A*X. If transpose = 'T', 
% then X is an (M x BLSZ) matrix and Y = A'*X.
% 
% OUTPUT OPTIONS:
% ---------------
%
% I.)   IRLBABLK(A) 
%       Displays the desired singular values in increasing order.   
%
% II.)  S = IRLBABLK(A) 
%       Returns the desired singular values in the vector S. 
%
% III.) [U,S,V] = IRLBABLK(A)
%       S is a diagonal matrix that contains the desired singular values in descending
%       order along the diagonal, the matrix V contains the corresponding "right" singular
%       vectors, and U contains the corresponding "left" singular vectors such that 
%       that A*V = U*S, A'*U = V*S, V'*V = I, and U'*U = I.  
%
% IV.)  [U,S,V,FLAG] = IRLBABLK(A)
%       This option returns the same as (III) plus a two dimensional array FLAG 
%       that reports if the algorithm converges and the number of matrix vector 
%       products. If FLAG(1) = 0 then this implies normal return: all desired 
%       singular values have converged. If FLAG(1) = 1 then the maximum number 
%       of iterations have been reached before all desired values converged. 
%       FLAG(2) contains the number of matrix vector products (A*X and A'*X is
%       recorded as 2*BLSZ matrix vector products.) used by the code. If the 
%       maximum number of iterations are reached (FLAG(1) = 1), then the matrices 
%       U, V, and S contain any singular pairs that have converged plus the last 
%       singular pair approximation for the pairs that have not converged.
%
% INPUT OPTIONS:
% --------------
%                                   
%       ... = IRLBABLK(A,OPTS)
%       OPTS is a structure containing input parameters. The input parameters can
%       be given in any order. The structure OPTS may contain some or all of the 
%       following input parameters. The string for the input parameters can contain
%       upper or lower case characters. 
%       
%  INPUT PARAMETER      DESCRIPTION    
%
%
%  OPTS.ADJUST       Initial number of vectors added to the K restart vectors. After
%                    vectors start to converge more vectors are added to help 
%                    increase convergence.
%                    DEFAULT VALUE    OPTS.ADJUST = 3
%
%  OPTS.AUG          Vectors used to augment the Krylov space. Choices are
%                    either Ritz vectors or harmonic Ritz vectors.
%                    DEFAULT VALUE    OPTS.AUG = 'HARM'  if SIGMA = 'SS'
%                                     OPTS.AUG = 'RITZ'  if SIGMA = 'LS'
%  
%  OPTS.BLSZ         Block size of the Lanczos bidiagonal matrix.        
%                    DEFAULT VALUE    BLSZ = 3               
%                        
%  OPTS.DISPS        Indicates if K approximate singular values and residuals are to be 
%                    displayed on each iteration. Set positive to display the values.
%                    DEFAULT VALUE   DISPS = 0 
%
%  OPTS.K            Number of desired singular values.             
%                    DEFAULT VALUE  K = 6
%
%  OPTS.MAXIT        Maximum number of iterations, i.e. maximum number of restarts.                           
%                    DEFAULT VALUE  MAXIT = 1000                      
%
%  OPTS.NBLS         Number of blocks in the Lanczos bidiagonal matrix. The program may increase
%                    NBLS to ensure certain requirements are satisfied. A warning message
%                    will be displayed if NBLS increases.                          
%                    DEFAULT VALUE    NBLS = 10
%
%  OPTS.REORTH       Three letter string specifying whether to use one-sided full reorthogonalization
%                    or two-sided. One-sided is performed only on the "short" vectors.
%                    Two-sided orthogonality is used when cond(A) estimated by cond(B) > 1/sqrt(eps).
%                    'ONE'  One-sided full reorthogonalization.
%                    'TWO'  Two-sided full reorthogonalization.
%                    DEFAULT VALUE  REORTH = 'ONE' 
%
%  OPTS.SIGMA        Two letter string specifying the location of the desired singular values.            
%                    'LS'  Largest singular values.
%                    'SS'  Smallest singular values.                 
%                    DEFAULT VALUE   SIGMA = 'LS'                                                         
%                                 
%  OPTS.TOL          Tolerance used for convergence. Convergence is determined when             
%                    || A*V - U*S || <= TOL*||A||. Where ||A|| is approximated by                                   
%                    the largest singular value of all projection matrices.  
%                    DEFAULT VALUE    TOL = 1d-6 
%                                                              
%  OPTS.V0           A matrix of approximate right singular vectors, if M >= N
%                    and sigma = 'LS' or 'SS'. If M < N and sigma = 'SS' then
%                    V0 must be a matrix of left singular vectors. This avoids
%                    computing zero eigenvalues of A'A. 
%                    DEFAULT VALUE  V0 = []
% 
%  DATE:  3/15/2006
%  VER:  1.0

%  AUTHORS:
%  James Baglama     University of Rhode Island  E-mail: jbaglama@math.uri.edu
%  Lothar Reichel    Kent State University       E-mail: reichel@math.kent.edu

% REFERENCES:
% 1.) J. Baglama and L. Reichel, "Augmented Implicitly Restarted Lanczos 
%     Bidiagonalization Methods", SIAM J. Sci. Comput., 27 (2005), pp. 19-42. 
% 2.) J. Baglama and L. Reichel, "Restarted Block Lanczos Bidiagonalization
%     Methods", submitted for publication.

% Incorrect number of output arguments requested.
if (nargout > 4 | nargout==2 ), error('ERROR: Incorrect number of output arguments.'); end

%----------------------------%
% BEGIN: PARSE INPUT VALUES. %
%----------------------------%

% No input arguments, return help.
if nargin == 0, help irlbablk, return, end

% Matrix A is stored in varargin{1}. Check type (numeric or character) and dimensions.
if (isstruct(varargin{1})), error('A must be a matrix.'), end
if ischar(varargin{1})
   if nargin == 1, error('Need dimension M for the (M x N) matrix A).'), end  
   m = varargin{2};
   if ~isnumeric(m) | length(m) ~= 1 
      error('Second argument M must be a numeric value.'); 
   end
   if nargin == 2, error('Need dimension N for the (M x N) matrix A).'); end  
   n = varargin{3}; 
   if ~isnumeric(n) | length(n) ~= 1 
      error('Third argument N must be a numeric value.'); 
   end
else
   [m,n] = size(varargin{1});
end

% Set all input options to default values.
adjust = 3; blsz=3; disps = 0;  K = 6; maxit = 1000; nbls = 10; sigma = 'LS';  
tol = 1d-6; aug = []; reorth='ONE'; V=[]; 

% Get input options from the structure.
if nargin > 1 + 2*ischar(varargin{1})
   options = varargin{2+2*ischar(varargin{1}):nargin};
   names = fieldnames(options);
   I = strmatch('ADJUST',upper(names),'exact');
   if ~isempty(I), adjust = getfield(options,names{I}); end
   I = strmatch('AUG',upper(names),'exact');
   if ~isempty(I), aug = upper(getfield(options,names{I})); end
   I = strmatch('BLSZ',upper(names),'exact');
   if ~isempty(I), blsz = upper(getfield(options,names{I})); end
   I = strmatch('DISPS',upper(names),'exact');
   if ~isempty(I), disps = getfield(options,names{I}); end
   I = strmatch('K',upper(names),'exact');
   if  ~isempty(I), K = getfield(options,names{I}); end
   I = strmatch('MAXIT',upper(names),'exact');
   if ~isempty(I), maxit = getfield(options,names{I}); end
   I = strmatch('NBLS',upper(names),'exact');
   if ~isempty(I), nbls = getfield(options,names{I}); end
   I = strmatch('REORTH',upper(names),'exact');
   if ~isempty(I), reorth = upper(getfield(options,names{I})); end
   I = strmatch('SIGMA',upper(names),'exact');
   if  ~isempty(I), sigma = upper(getfield(options,names{I})); end
   I = strmatch('TOL',upper(names),'exact');
   if ~isempty(I), tol = getfield(options,names{I}); end
   I = strmatch('V0',upper(names),'exact');
   if ~isempty(I), V = getfield(options,names{I}); end
end 

% Check type of input values and output an error message if needed.
if (~isnumeric(adjust) | ~isnumeric(disps)   | ~isnumeric(K)  | ~isnumeric(maxit) | ...
    ~isnumeric(nbls)    | ~isnumeric(tol) | ~ischar(sigma) | ~ischar(reorth) | ~isnumeric(blsz))
   error('Incorrect type for input value(s) in the structure.');
end

% Check the values of sigma.
if length(sigma)  ~= 2, error('SIGMA must be LS or SS'); end
if (~strcmp(sigma,'SS') & ~strcmp(sigma,'LS') )
   error('SIGMA must be a two letter string LS or SS.'); 
end

% Check the values of reorth.
if length(reorth) ~=3, error('REORTH must be the three letter string ONE or TWO.'); end
if (~strcmp(reorth,'ONE') & ~strcmp(reorth,'TWO') )
   error('REORTH must be the three letter string ONE or TWO.'); 
end

% Interchange M and N so that size(A'A) = min(m,n). This avoids
% finding zero values when searching for the smallest singular
% values. No interchange needed if SIGMA = 'LS'.
interchange = 0; if n > m & strcmp(sigma,'SS'), t=m; m=n; n=t; interchange = 1; end 

% Increase the number of desired values by adjust to help increase convergence. K
% is re-adjusted as vectors converge. This is only an initial value of K.
K_org = K; K = K + adjust; 

% Check for input errors in the data structure.
if K       <= 0,    error('K must be a positive value.'),    end
if blsz    <= 0,    error('BLSZ must be a positive value.'), end
if K  >  min(n,m),  error('K must be less than min(n,m) + %g.',adjust),  end   
if nbls     <= 1,   error('NBLS must be greater than 0.'),   end
if tol     <  0,    error('TOL must be non-negative.'),      end
if maxit   <= 0,    error('MAXIT must be positive.'),        end
if blsz*nbls >= min(n,m)
   nbls = floor((min(n,m)-blsz - 0.1)/blsz);
   warning(['Changing NBLS to ',num2str(nbls)]);
   if blsz*nbls - K -  blsz  < 0
      adjust = 0; K = blsz*nbls - blsz;
      warning(['Changing adjust to ',num2str(adjust)]);
      warning(['Changing K to ',num2str(K)]);
   end
end
if blsz*nbls - K - blsz  < 0
   nbls = ceil((K+blsz)/blsz+blsz+0.1);
   warning(['Changing NBLS to ',num2str(nbls)]);
end
if blsz*nbls >= min(n,m)
   nbls = floor((min(n,m)-blsz - 0.1)/blsz);
   adjust = 0; K = blsz*nbls - blsz; 
   warning(['Changing adjust to ',num2str(adjust)]);
   warning(['Changing K to ',num2str(K)]);
   warning(['Changing NBLS to ',num2str(nbls)]);
end

% Preallocate memory for W and F. These matrices are full and resizing will cause
% an increase in cpu time.
W = zeros(m,blsz*nbls); F = zeros(n,blsz);

% If starting matrix V0 is not given then set starting matrix V0 to be a 
% (n x 1) matrix of normally distributed random numbers.
if isempty(V)
  V = zeros(n,blsz*nbls); % Preallocate memory for V. This matrix is a full matrix.
  V(:,1:blsz) = randn(n,blsz); 
else 
  if (size(V,2) ~= blsz), error('Incorrect size of starting matrix V0.'), end
  V(:,blsz+1:blsz*nbls) = zeros(n,blsz*nbls-blsz);
end
if ~isnumeric(V), error('Incorrect starting matrix V0.'), end
if (size(V,1) ~= n), error('Incorrect size of starting matrix V0.'), end 

% Determine which vectors to use for augmenting.
if isempty(aug)
  if strcmp(sigma,'LS'), aug = 'RITZ'; else aug = 'HARM'; end
else
  if length(aug) ~= 4, error('Unknown value for AUG. AUG must be the four letter string RITZ or HARM'); end
  if ~ischar(aug), error('Unknown value for AUG. AUG must be the four letter string RITZ or HARM'); end
  if (~strcmp(aug,'RITZ') & ~strcmp(aug,'HARM') )
     error('Unknown value for AUG. AUG must be the four letter string RITZ or HARM'); 
  end
end

% Set tolerance to machine precision if tol < eps.
if tol < eps, tol = eps; end

%--------------------------%
% END: PARSE INPUT VALUES. %
%--------------------------%

%-----------------------------------------------------------%
% BEGIN: DESCRIPTION AND INITIALIZATION OF LOCAL VARIABLES. %
%-----------------------------------------------------------%

% Initialization and description of local variables.   
B = [];              % Block bidiagonal matrix.
Bsz =[];             % Size of the block bidiagonal matrix.
conv = 'F';          % Boolean to determine if all desired singular values have converged.
eps23 = eps^(2/3);   % Two thirds of eps, used for Smax. Avoids using zero.
I=[];                % Used for indexing.
iter = 1;            % Main loop iteration count.
J=[];                % Used for indexing.
mprod = 0;           % The number of matrix vector products with A and A'.
R_F =[];             % The R matrix from the QR factorization of the residual matrix F.
sqrteps = sqrt(eps); % Square root of machine tolerance used in convergence testing.
Smax = 1;            % Holds the maximum value of all computed singular values of B est. ||A||_2.
Smin=[];             % Holds the minimum value of all computed singular values of B est. cond(A).
SVTol = min(sqrteps,tol);  % Tolerance to determine when a singular vector has converged.
S_B =[];             % Singular values of B.
U_B =[];             % Left singular vectors of B.
V_B =[];             % Right singular vectors of B.
V_B_last=[];         % Holds the last row of the modified V_B.
S_B2 =[];            % Singular values of [B (E_m*R_F')].
U_B2 =[];            % Left singular vectors of [B (E_m*R_F')].
V_B2 =[];            % Right singular vectors of [B (E_m*R_F')].

%--------------------------------------------------------------------%
% END: DESCRIPTION AND INITIALIZATION OF LOCAL AND GLOBAL VARIABLES. %
%--------------------------------------------------------------------%

%----------------------------%
% BEGIN: MAIN ITERATION LOOP %
%----------------------------%

while (iter <= maxit)
    
    % Compute the Lanczos bidiagonalization decomposition.
    [V,W,F,B,mprod] = ablanzbd(varargin{1},V,W,F,B,K,interchange,nbls,blsz,n,...
                               m,mprod,SVTol*Smax,reorth,iter);

    % Determine the size of the bidiagonal matrix B.
    Bsz = size(B,1);  

    % Compute the QR-factorization of F.
    [F,R_F] = qr(F,0);

    % Compute singular triplets of B. MATLAB's svds orders the singular values
    % largest to smallest.
    [U_B,S_B,V_B] = svd(B); S_B = diag(S_B);
    
    % Estimate ||A|| using the largest singular value over all iterations
    % and estimate the cond(A) using approximations to the largest and smallest 
    % singular values. If a small singular value is less than sqrteps use only Ritz
    % vectors to augment and require two-sided reorthogonalization, i.e. set REORTH = 'TWO'.
    if iter==1
       Smax = S_B(1); Smin = S_B(Bsz);
    else 
       Smax = max(Smax,S_B(1)); Smin = min(Smin,S_B(Bsz));
    end
    Smax = max(eps23,Smax); % Automatic switch from aug = 'HARM' to aug = 'RITZ'.
    if Smin/Smax < sqrteps, reorth = 'TWO'; aug = 'RITZ'; 'cross', iter, end
   
    % Re-order the singular values accordingly. MatLab's SVD orders the 
    % singular values largest to smallest.
    if strcmp(sigma,'SS'), I = Bsz:-1:1; U_B = U_B(:,I); V_B = V_B(:,I); S_B = S_B(I); end   
   
    % Compute the residuals for the singular values.
    R = R_F*U_B(Bsz-(blsz-1):Bsz,:); residuals = (sum(R.^2, 1)).^(0.5); 
    
    % Convergence tests and displays.
    [conv,U_B,S_B,V_B,K] = convtests(Bsz,blsz,disps,tol,K_org,U_B,S_B,V_B,residuals,iter,K,SVTol,Smax);

    % If all desired singular values converged then exit main loop.
    if strcmp(conv,'T'),  break,  end;

    % Reached maximum number of iterations, exit main loop.
    if iter >= maxit, break, end;
    
    % Compute the starting vectors and first block B(1:K,1:K+1). 
    if strcmp(aug,'HARM')

       % Update the SVD of B to be the SVD of [B (E_m*R_F')].
       [U_B2,S_B,V_B2] = svd([diag(S_B) R']); S_B = diag(S_B);   
       if strcmp(sigma,'SS')
          U_B2 = U_B2(:,I); V_B2 = V_B2(:,I); S_B = S_B(I);
       end
       U_B = U_B*U_B2;
       EB_blsz = rot90(flipud(eye(Bsz+blsz,blsz)))'; 
       V_B = [[V_B; zeros(blsz,Bsz)] EB_blsz]*V_B2;
       V_B_last = V_B(Bsz+1:Bsz+blsz,1:K); % Set equal to the last BLSZ rows of V_B.
       EB = rot90(flipud(eye(Bsz,blsz)))';
       s = B\EB*R_F'; 
       V_B = V_B(1:Bsz,:) + s*V_B(Bsz+1:Bsz+blsz,:);
     
       % Vectors V_B are not orthogonal.
       [V_B,R] = qr([ [V_B(:,1:K); zeros(blsz,K)]  [-s; eye(blsz)] ],0);
       V(:,1:K+blsz) = [V(:,1:Bsz) F]*V_B; 
  
       % Update and compute the K x K+BLSZ part of B.
       B = diag(S_B(1:K))*triu((R(1:K+blsz,1:K) + R(:,K+1:K+blsz)*V_B_last)');       

    else
       V(:,1:K+blsz) = [V(:,1:Bsz)*V_B(:,1:K) F];
       B = [diag(S_B(1:K)), R(:,1:K)'];
    end

    % Compute left approximate singular vectors.
    W(:,1:K) = W(:,1:Bsz)*U_B(:,1:K);

    % Update the main iteration loop count.
    iter = iter + 1;

end 

%--------------------------%
% END: MAIN ITERATION LOOP %
%--------------------------%

%-----------------------%
% BEGIN: OUTPUT RESULTS %
%-----------------------%

% Test to see if maximum number of iterations have been reached.
FLAG(1) = 0; if iter >= maxit FLAG(1) = 1; end; 

% Output option I: Display singular values only.
if (nargout == 0)
  if FLAG(1)==1, disp('Maximum number of iterations exceeded.'); end
  SingularValues = S_B(1:K_org)
end

% Output option II: Set singular values equal to output vector.  
if (nargout == 1), varargout{1} = S_B; end

% Output option III and IV: Output singular triplets (U,S,V), where
% A*V = U*S and A'*U = V*S.
if nargout > 1
   if interchange
      varargout{1} = V(:,1:Bsz)*V_B(:,1:K_org); 
      varargout{2} = diag(S_B(1:K_org)); 
      varargout{3} = W(:,1:Bsz)*U_B(:,1:K_org);
   else
      varargout{1} = W(:,1:Bsz)*U_B(:,1:K_org); 
      varargout{2} = diag(S_B(1:K_org)); 
      varargout{3} = V(:,1:Bsz)*V_B(:,1:K_org);
   end
   if nargout == 4, FLAG(2) = mprod; varargout{4} = FLAG; end
end

%---------------------%
% END: OUTPUT RESULTS %
%---------------------%

%------------------------------------------------%
% BEGIN: LANCZOS BIDIAGONALIZATION DECOMPOSITION %
%------------------------------------------------%

function [V,W,F,B,mprod] = ablanzbd(A,V,W,F,B,K,interchange,nbls,blsz,n,...
                                    m,mprod,SVTol,reorth,iter)
% Computes the block Lanczos bidiagonalization decomposition
%  
%  A*V  = W*B
%  A'*W = V*B' + F*E^T
%  
% with full reorthogonalization. The matrix A can be passed as a numeric 
% matrix or as a filename. If the M x N matrix A is a filename then
% Y = Afunc(X,M,N,BLSZ,'transpose'). If transpose = 'F', then X is an (N x BLSZ)
% matrix and Y = A*X. If transpose = 'T', then X is an (M x BLSZ) matrix
% and Y = A'*X.

% James Baglama
% DATE: 3/15/2006

% Initialization of main loop count J.
J = 1:blsz; 

% First orthogonalization of starting vectors via QR factorization.
if iter == 1
   [V,R] = qr(V,0); B=[];
else
   J = J + K;
end

% Matrix A product with vector(s) V, (W=A*V).
if interchange
   if ischar(A)
     W(:,J) = feval(A,V(:,J),m,n,blsz,'T');
  else
     W(:,J) = (V(:,J)'*A)';
  end
else
  if ischar(A)
     W(:,J) = feval(A,V(:,J),m,n,blsz,'F');
  else
     W(:,J) = A*V(:,J);
  end
end

% Count the number of matrix vector products.
mprod = mprod + blsz;

% Input vectors are singular vectors and AV(:,J) which must be orthogonalized.
if iter ~= 1
   W(:,J) = orthog(W(:,J),W(:,1:J(blsz)-blsz));
end

% Compute the QR factorization of W.
[W(:,J),S,E] = qr(W(:,J),0); 

% Check for linearly dependent vectors.
stol = max(SVTol,eps*abs(S(1,1)));
if abs(S(1,1)) > stol
   rat_diag = abs(S(blsz,blsz)/S(1,1));
else
   rat_diag = 0; % All vectors in the block are zero.
end 
EE(E)=[1:blsz];  S = S(:,EE);
if rat_diag <= stol
   W(:,J) = orthog(W(:,J),W(:,1:J(blsz)-blsz));
   [W(:,J),S_2] = qr(W(:,J),0); S = S_2*S; 
end

% Begin of main iteration loop for the block Lanczos bidiagonalization decomposition.
while (J(blsz) <= blsz*nbls)

   % Matrix A' product with vector(s), (F = A'*W). 
   if interchange
      if ischar(A)
         F = feval(A,W(:,J),m,n,blsz,'F');
      else
         F = A*W(:,J);  
      end
   else
      if ischar(A)
        F = feval(A,W(:,J),m,n,blsz,'T');
      else
        F = (W(:,J)'*A)';
      end
   end
     
   % Count the number of matrix vector products.
   mprod = mprod + blsz;
   
   % One step of the block classical Gram-Schmidt process.
   F = F - V(:,J)*S';

   % Full Reorthogonalization step. "Short vectors".
   F = orthog(F,V(:,1:J(blsz)));

   if J(blsz)+blsz <= blsz*nbls

       % Compute the QR factorization of F.
       [V(:,J+blsz),R,E] = qr(F,0);
      
       % Check for linearly dependent vectors.
       stol = max(SVTol,eps*abs(R(1,1)));
       if abs(R(1,1)) > stol
          rat_diag = abs(R(blsz,blsz)/R(1,1));
       else
           rat_diag = 0; % All vectors in the block are numerically zero.
       end 
       EE(E)=[1:blsz];  R = R(:,EE);
       if rat_diag <= stol
         V(:,J+blsz) = orthog(V(:,J+blsz),V(:,1:J(blsz)));
         [V(:,J+blsz),R_2] = qr(V(:,J+blsz),0); R = R_2*R;
       end 
           
       % Compute block bidiagonal matrix B. 
       if isempty(B) 
         B = [S R'];
       else
         B = [B zeros(J(1)-1,blsz); zeros(blsz,J(1)-1) S R'];
       end 
    
       % Matrix A product with vector(s), (W=A*V).
       if interchange
         if ischar(A)
           W(:,J+blsz) = feval(A,V(:,J+blsz),m,n,blsz,'T');
         else
           W(:,J+blsz) = (V(:,J+blsz)'*A)';
         end
       else
         if ischar(A)
            W(:,J+blsz) = feval(A,V(:,J+blsz),m,n,blsz,'F');
         else
            W(:,J+blsz) = A*V(:,J+blsz);
         end
       end

       % Count the number of matrix vector products.
       mprod = mprod + blsz;

       % One step of the block classical Gram-Schmidt process.
       W(:,J+blsz) =  W(:,J+blsz) - W(:,J)*R';

      % Full Reorthogonalization step. "Long vectors"
      if ( iter == 1 | strcmp(reorth,'TWO') )
         W(:,J+blsz) = orthog(W(:,J+blsz),W(:,1:J(blsz)));
      end
   
      % Compute the QR factorization of W.
      [W(:,J+blsz),S,E] = qr(W(:,J+blsz),0);
 
      % Check for linearly dependent vectors.
      stol = max(SVTol,eps*abs(S(1,1)));
      if abs(S(1,1)) > stol
        rat_diag = abs(S(blsz,blsz)/S(1,1));
      else
        rat_diag = 0; % All vectors in the block are zero.
      end 
      EE(E)=[1:blsz];  S = S(:,EE);
      if rat_diag <= stol 
         W(:,J+blsz) = orthog(W(:,J+blsz),W(:,1:J(blsz)));
        [W(:,J+blsz),S_2] = qr(W(:,J+blsz),0); S = S_2*S;
      end      

   else
    
    % Add last block to matrix B
     B = [B; zeros(blsz,J(1)-1) S]; 

   end

    % Update iteration count.
    J = J + blsz;

end

function Y = orthog(Y,X)
% Orthogonalize vectors Y against vectors X.
dotY = (Y'*X)'; Y = Y - X*dotY; 

 
%----------------------------------------------------%
% END: BLOCK LANCZOS BIDIAGONALIZATION DECOMPOSITION %
%----------------------------------------------------%

%--------------------------%
% BEGIN: CONVERGENCE TESTS %
%--------------------------%

function [conv,U_B,S_B,V_B,K] = convtests(Bsz,blsz,disps,tol,K_org,U_B,S_B,V_B,residuals,iter,K,SVTol,Smax)
% This function checks the convergence of singular values. 

% James Baglama
% DATE: 3/15/2006

% Set convergence to false.
conv = 'F';  

% Output intermediate results.
 if disps ~= 0 
   disp(sprintf('  Singular Value     Residual       Iteration: %d',iter));
   S = sprintf('%15.5e %15.5e \n',[S_B(1:K_org)';residuals(1:K_org)]);
   disp(S); disp(' '); disp(' ');
 end

% Check for convergence.
Len_res = length(find(residuals(1:K_org)' < tol*Smax));
if Len_res == K_org 
   U_B = U_B(:,1:K_org); V_B = V_B(:,1:K_org); S_B = S_B(1:K_org);
   
   % Set convergence to true.
   conv = 'T'; return;
   
end

% Adjust K to include more vectors as the number of vectors converge.
Len_res = length(find(residuals(1:K_org)' < SVTol*Smax));
K = max(K, K_org + Len_res); if K > Bsz - blsz - 3, K = Bsz - blsz - 3; end

%------------------------%
% END: CONVERGENCE TESTS %
%------------------------%