function [x,flag,Arrnorm] = alsqr(A,b,x,OPTS)
                                             
% ALSQR is an augmented preconditioned LSQR method. The 
% preconditioner is computed via an augmented Lanczos 
% bidiagonalization method with harmonic Ritz vectors. The program then
% runs LSQR on the preconditioned system. The ALSQR method solves the 
% least squares problem: 
%
%                  min || b - A*x || 
%                   x
%
% where A is an (m x n) matrix, b an (m x 1) vector, and x is an (n x 1) 
% vector. b is not assumed to be in span(Col(A)). ||  || always denotes
% 2 norm.
%
%-------------------------------------------------------------------------
%
% [x,FLAG,Arrnorm] = ALSQR(A,b,x)
%
% The first input argument into ... = alsqr(A,b,x) can be a numeric 
% matrix A or an M-file 'Afunc'. If the m x n matrix A is a filename then 
% z = Afunc(x,1) computes z = A*x and z = Afunc(x,2) computes z = A'*x.
% The second input is the vector b and third input is the vector x, an
% initial approximation to min || b - A*x ||, typically chosen to be 0.
%
%-------------------------------------------------------------------------
%
% OUTPUT OPTIONS:
% ---------------
%    
% I.)   x = ALSQR(A,b,x) 
%       Returns only an approximate solution vector x of min || b - A*x ||.                                              
%
% II.)  [x,FLAG,Arrnorm] = ALSQR(A,b,x,OPTS)
%       x -      An approximate solution of min || b - A*x ||.
%
%       FLAG  - Four dimensional array 
%               (r = b-A*x and r_0 = b-A*x_0 where x is the current 
%                approximation to min || b - A*x || and x_0 is the initial 
%                approximation on input.):
%               FLAG(1) = 1  Stopping rule is satisfied:   
%                            ||A'*r|| <= LTOL*||A'*r_0||
%                       = 0  Stopping rule is not satisfied.
%
%               FLAG(2) = The approximate condition number, cond(A), 
%                         of the matrix A.
%
%               FLAG(3) = Number of iterations of part I building the 
%                         preconditioner.
%
%               FLAG(4) = Number of iterations of part II nonrestarted LSQR.
%
%      Arrnorm - Three dimensional array such that 
%                Arrnorm(j,1) = ||A'*(b - A*x)|| = ||A'*r||
%                Arrnorm(j,2) = ||b - A*x|| = ||r||
%                and Arrnorm(j,3) is the total number of matrix-vector 
%                products with A and A' required to compute the solution x(j).
%
%-------------------------------------------------------------------------
%
% INPUT OPTIONS:
% --------------
%                                   
%       [x,FLAG,Arrnorm] = ALSQR(A,b,x,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.
%       alsqr(A,b) without the structure OPTS will use all of the default
%       values for OPTS.
%
%  INPUT PARAMETER      DESCRIPTION    
%
%
%  OPTS.ADJUST       Number of vectors to add to the k augmenting vectors.
%                    DEFAULT VALUE   OPTS.ADJUST = 40 
%
%  OPTS.HTOL         Tolerance used for convergence of harmonic Ritz vector
%                    approximations to the k smallest singular triplets in
%                    preconditioning stage, part I. Convergence is determined              
%                    sqrt(|| A'*U - V*S ||^2 + || A*V - U*S ||^2) <= HTOL*||A||
%                    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. 
%                    DEFAULT VALUE    HTOL = 1d-3
%
%  OPTS.K            Number of harmonic Ritz vectors to augment Krylov 
%                    subspace to precondition the system. Setting k = 0
%                    will call the nonrestarted LSQR method only.
%                    DEFAULT VALUE   K = 20
%
%  OPTS.LTOL         Tolerance used for stopping rule.
%                    ||A'*r|| <= LTOL*||A'*r_0||
%                    DEFAULT VALUE    LTOL = 1d-6 
%
%  OPTS.MAXITP       Maximum number of iterations in building the 
%                    preconditoner. Part I.                   
%                    DEFAULT VALUE   MAXITP = 1000 
%
%  OPTS.MAXITL       Maximum number of iterations of nonrestarted LSQR.
%                    Part II.
%                    DEFAULT VALUE   MAXITL = min(m,n). 
%
%  OPTS.M_B          Size of the Lanczos bidiagonal matrix. The program may 
%                    increase M_B to ensure certain requirements are 
%                    satisfied. A warning message will be displayed if M_B 
%                    increased.                          
%                    DEFAULT VALUE   M_B = 120
%
%  OPTS.M_REORTH     Number of vectors to reorthogonalize back against 
%                    in part II, nonrestarted LSQR stage. Set M_REORTH to
%                    MAXITL performs reorthogonalization of all vectors 
%                    computed. Does not apply Part I. M_REORTH must be set
%                    to be >= M_B.
%                    DEFAULT VALUE    M_REORTH = M_B
%
%  OPTS.REORTH012    Part I preconditioning stage: Full reorthogonalization 
%                    (one or two sided) is always performed on the M_B vectors 
%                    during the preconditioning stage. Part II nonrestarted 
%                    LSQR stage: Either no reorthogonalization, one-sided 
%                    or two-sided reorthogonalization is performed. 
%                    REORTH012 = 0 One-sided reorthogonalization 
%                                  preconditioning stage  (Part I) and no 
%                                  reorthogonalization for nonrestarted 
%                                  LSQR stage (Part II).
%                    REORTH012 = 1 One-sided reorthogonalization 
%                                  preconditioning stage  (Part I) and one-sided 
%                                  reorthogonalization for nonrestarted 
%                                  LSQR stage (Part II).
%                    REORTH012 = 2 Two-sided reorthogonalization 
%                                  preconditioning stage  (Part I) and two-sided 
%                                  reorthogonalization for nonrestarted 
%                                  LSQR stage (Part II).
%                    REORTH012 = 3 Two-sided reorthogonalization 
%                                  preconditioning stage  (Part I) and one-sided 
%                                  reorthogonalization for nonrestarted 
%                                  LSQR stage (Part II).  
%                    DEFAULT VALUE REORTH012 = 0 
%
%-------------------------------------------------------------------------
%
%  DATE:  4/20/12
%  VER:  1.0
%
%  AUTHORS:
%  James Baglama     
%  University of Rhode Island 
%  EMAIL: jbaglama@math.uri.edu
% 
%  Lothar Reichel
%  Kent State University
%  EMAIL: reichel@math.kent.edu
%
%  Dan Richmond
%  University of Rhode Island 
%  EMAIL: dan@math.uri.edu 
%
% REFERENCES:
% 1.) J. Baglama, L.Reichel, and D. Richmond, "An Adaptively Preconditioned LSQR Algorithm".
% 2.) J. Baglama and L. Reichel, "Augmented Implicitly Restarted Lanczos 
%     Bidiagonalization Methods", SIAM J. Sci. Comput., 27 (2005), pp. 19-42. 
%
%-------------------------------------------------------------------------

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

% Wrong number of input or output values.
if (nargin < 3 || nargin > 4 || nargout < 1 || nargout==2 || nargout >3 )
    error('alsqr: Incorrect number of input or output arguments.'); 
end

% Matrix A can be input as a numeric matrix or as a function handle.
if isa(A,'numeric')
    [m,n] = size(A); Anumeric = true;
elseif isa(A,'function_handle')
    m = size(b,1); n = size(x,1); Anumeric = false;
else
   error('alsqr: Unknown type for matrix A.');
end

% Set all input options to default values.
adjust = 40; k = 20; htol = 1d-3; ltol = 1d-6; maxitp = 1000; 
maxitl=min(m,n); m_b = 120; m_reorth = m_b; reorth012 = 0; 

% Get user input options from the structure.
if nargin == 4
   if ~isa(OPTS,'struct')
            error('alsqr:OptionsNotStructure',...
               'Options argument must be a structure.')
   end
   names = fieldnames(OPTS);
   I = strcmpi('ADJUST',names); I=find(I>0);
   if ~isempty(I), adjust = OPTS.(names{I}); end
   I = strcmpi('K',names);I=find(I>0);
   if  ~isempty(I), k = OPTS.(names{I}); end
   I = strcmpi('HTOL',names); I=find(I>0);
   if ~isempty(I), htol = OPTS.(names{I}); end
   I = strcmpi('LTOL',names); I=find(I>0);
   if ~isempty(I), ltol = OPTS.(names{I}); end
   I = strcmpi('MAXITP',names); I=find(I>0);
   if ~isempty(I), maxitp = OPTS.(names{I}); end
   I = strcmpi('MAXITL',names); I=find(I>0);
   if ~isempty(I), maxitl = OPTS.(names{I}); end
   I = strcmpi('M_B',names); I=find(I>0);
   if ~isempty(I), m_b = OPTS.(names{I}); m_reorth = m_b; end
   I = strcmpi('M_REORTH',names); I=find(I>0);
   if ~isempty(I), m_reorth = OPTS.(names{I}); end
   if m_reorth < m_b, error('alsqr: Requires m_reorth >= m_b'); end
   I = strcmpi('REORTH012',names); I=find(I>0);
   if ~isempty(I), reorth012 = OPTS.(names{I}); end
end 

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

%-------------------------------%
% BEGIN: SETUP INPUT VARIABLES. %
%-------------------------------%
    
% Preallocate memory for all large matrices and vectors. 
% Resizing will cause an increase in cpu time.
if m_reorth <= m_b+1
  Q = zeros(m,m_b+1);    % Full matrix, left singular (harmonic) vectors of A.
  P = zeros(n,m_b+1);    % Full matrix, right singular vectors of A.
else
  Q = zeros(m,m_reorth); P = zeros(n,m_reorth);
end  
w = zeros(n,1);        % Storage vector required for LSQR.
  
% Setup the ouput variable Arrnorm.
Arrnorm = [ones(maxitp+maxitl,2) zeros(maxitp+maxitl,1)];

% Compute the initial residual if x is not zero.
% Q(:,1) is the starting vector for the routine.
% The matrix vector computed here is not counted
% in Arrnorm.
if nnz(x) ~= 0
    if Anumeric
        Q(:,1) = b - A*x;    
    else
        Q(:,1) = b - A(x,1); 
    end
else
   Q(:,1) = b;
end   

% Initialize the output value flag.
flag = zeros(1,4);

%-----------------------------%
% END: SETUP INPUT VARIABLES. %
%-----------------------------%


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

% Initialization and description of local variables.
alpha=[];              % alpha_{m+1} in the Lanczos decomposition.
B = [];                % Bidiagonal matrix.    
hcount=0;              % Count of converged harmonic Ritz vectors.
k_org = k;             % Holds initial value of k.
nrcount = 0;           % Count for residual norms. 
resHarm = [];          % Residual norms of harmonic Ritz vectors.
piter = 1;             % Main loop iteration count for building preconditioner Part I.
rnorm = [];            % Holds ||B-A*X(J))|| and matrix-vector products.
f = [];                % Initial right hand side. Small vector, maximum size m_b.
f_last = [];           % Last element of the right hand side, used to compute residual.
s = [];                % Place holder.
Smax = -1;             % Holds the maximum value of all computed singular values of B.
Smin = Inf;            % Holds the minimum value of all computed singular values of B.
S_B = [];              % Singular values of B.
U_B = [];              % Left singular vectors of B.
U_old = [];            % Place holder. Left singular values of B, for updating harmonic Ritz vectors. 
U_B_last = [];         % Holds the last column of U_B.
V_B = [];              % Right singular vectors of B.
%----------------------------------------------------%
% END: DESCRIPTION AND INITIALIZATION OF  VARIABLES. %
%----------------------------------------------------%

%--------------------------------------------------------------%
% BEGIN: PART I BUILDING PRECONDITIONER AND UPDATING SOLUTION  %
%--------------------------------------------------------------%

% Increase the number of desired values by adjust to help increase convergence. K
% is re-adjusted as vectors converge. If k = 0, then no augmenting 
% and the nonrestarted LSQR is called.
if k > 0, k = k + adjust; end 
 

% Initial call to function augbidiaglsqr. Set k = 0 (no vectors for augmenting).
% On exit, if orgional k > 0 we will have B is a m_b+1 x m_b, Q is 
% a m x m_b+1, P is a n x m_b+1.
if k > 0
  [alpha,Arrnorm,B,flag,f_last,nrcount,P,Q,w,x] = augbidiaglsqr(A,...
   Anumeric,Arrnorm,B,f,flag,0,ltol,m_b,m_b,m,n,nrcount,P,Q,reorth012,w,x);
end   
  
% Main iteration loop to build up the preconditioner.  
while (piter <= maxitp && ~flag(1) && k > 0)         
      
     % Compute SVD of the rectangular lower bidiagonal matrix B.
     % Matlab orders largest to smallest. 
     % B is a m_b+1 x m_b, U_B is a m_b+1 x m_b+1, V_B is a m_b x m_b
     [U_B,S_B,V_B] = svd(B);
     U_B_last = U_B(:,m_b+1); % Last column of U_B
     U_B      = U_B(:,1:m_b); S_B=S_B(1:m_b,:); %  Reset to economy size or else reordering fails.
     U_B      = U_B(:,end:-1:1); % Columns in reverse order.
     S_B(:)   = S_B(end:-1:1);   % Diagonal with smallest singular value at S_B(1,1). 
     V_B      = V_B(:,end:-1:1); % Writing the columns in reverse order.

     % Estimate ||A|| and condition number of A using the largest and smallest singular 
     % values over all iterations.  
     Smax = max(Smax,S_B(m_b,m_b)); Smin = min(Smin,S_B(1,1));  
     flag(2) = Smax/Smin; % Approximate condition number of A. 
      
     % Holds the old value of U_B that is used to compute new matrix B.
     U_old = U_B; 
    
     % Compute the harmonic Ritz vectors.  
     s   = -U_B_last(1:m_b,1)/(U_B_last(m_b+1,1)); 
     U_B =  U_B(1:m_b,:)+s*U_B(m_b+1,:);          
      
     % Compute QR of U_B and the residual vector using MATLAB's internal qr function.
     [U_B,~] = qr([ [U_B(:,1:k); zeros(1,k)] [-s; 1] ],0);

     % Compute the vector f.
     f = f_last*[-s; 1];
      
     % Update U, V and f.
     Q(:,1:k+1) = Q(:,1:m_b+1)*U_B;
     P(:,1:k+1) = [P(:,1:m_b)*V_B(:,1:k) P(:,m_b+1)];
     f = U_B'*f;
      
     hcount=0; % Reset the counter of converged harmonic Ritz vectors. 
     for j=1:k_org
        if (sqrt( (S_B(j,j)^2)*(norm(U_old(:,j)-U_B(:,j)))^2  + ...
                  (norm(B'*U_B(:,j)-S_B(j,j)*V_B(:,j)))^2     + ...
                  (alpha*U_B(m_b+1,j))^2) < htol*Smax)   
                     hcount = hcount + 1; 
        end
     end
      
     % Update B and alpha. Matlab command TRIL is used to ensure
     % triangular part of B.
     B = tril([S_B(1:k,1:k)*U_old(:,1:k)'*U_B; alpha*U_B(m_b+1,:)]');
      
     % Break from preconditioning stage part I and continue to nonrestarted
     % lsqr part II if all harmonic Ritz vectors have converged
     if (hcount >= k_org || piter >= maxitp), break; end
      
     % Update iteration count.
     piter = piter + 1; 
     flag(3) = piter;
      
     % Call the restarted augmented harmonic method to build the space,
     % approximate the solution, and update the residual norms.
     % On exit we will have B is a m_b+1 x m_b, Q is a m x m_b+1, 
     % P is a n x m_b+1 
     [alpha,Arrnorm,B,flag,f_last,nrcount,P,Q,w,x] = augbidiaglsqr(A,...
      Anumeric,Arrnorm,B,f,flag,k,ltol,m_b,m_b,m,n,nrcount,P,Q,reorth012,w,x);
         
end      
 
%------------------------------------------------------------%
% END: PART I BUILDING PRECONDITIONER AND UPDATING SOLUTION  %
%------------------------------------------------------------%

%-----------------------------------%
% BEGIN: PART II NONRESTARTED LSQR  %
%-----------------------------------%

% Set the number of nonrestarted LSQR to zero.
flag(4) = 0;

% Reset reorth to be one-sided for nonrestarted LSQR.
if reorth012 == 3, reorth012 = 1; end 

% Call the nonrestarted lsqr part of function augbidiaglsqr.
if ~flag(1)
  [alpha,Arrnorm,B,flag,f_last,nrcount,P,Q,w,x] = augbidiaglsqr(A,...
   Anumeric,Arrnorm,B,f,flag,k,ltol,maxitl,m_reorth,m,n,nrcount,P,Q,reorth012,w,x);
end

%---------------------------------%
% END: PART II NONRESTARTED LSQR  %
%---------------------------------%


%-------------------------%       
% FUNCTION: AUGBIDIAGLSQR %   
%-------------------------%

function [alpha,Arrnorm,B,flag,f_last,nrcount,P,Q,w,x] = augbidiaglsqr(A,...
          Anumeric,Arrnorm,B,f,flag,k,ltol,maxitl,m_b,m,n,nrcount,P,Q,reorth012,w,x) 
% 
% Function computes the Lanczos Bidiagonalization and updates LSQR solution
% and residual. The function also can continue onto nonrestarted LSQR. 
% INPUT:
%       A - (m x n) matrix
%       Anumeric - logic variable indicates whether A is a function call or
%                  numeric variable.
%       Arrnorm  - (3 x 1) array contains the ||A'*r||, ||r|| and # 
%                  of matrix-vector products.
%       B  - (k+1 x k) bidiagonal matrix used in the Lanczos bidiagonalization
%       f  - (k+1 x 1) vector right hand side used for least squares
%            solution.
%       k  - scalar number of augmenting vectors.
%     ltol - scalar tolerance used for convergence of LSQR
%   maxitl - scalar maximum number of iteration for LSQR, set to m_b during
%            the precondiioning stage.
%      m_b - scalar size of the lanczos bidiagonalization and the number of
%            reorthogonalization vectors during the LSQR stage.
%       m  - scalar number of rows of A
%       n  - scalar number of columns of A
%  nrcount - scalar count used for Arrnorm
%        P - (n x k+1) orthogonal matrix used in the Lanczos bidiagonalization
%        Q - (m x k+1) orthogonal matrix used in the Lanczos bidiagonalization
% reorth012- scalar indiactes no, one or two sided reorthogonalization
%       w  - (n x 1) vector used to update soliution x in LSQR
%       x  - (n x 1) vector the LSQR solution
%
% OUTPUT:
%    alpha - scalar holds the last diagonal element of B.
% Arrnorm  - (3 x 1) updated array that contains the ||A'*r||, ||r|| 
%                   and # of matrix-vector products. 
%       B  - (m_b+1 x m_b) bidiagonal matrix used in the Lanczos bidiagonalization
%     flag - (4 x 1) array holds convegence indicator
%   f_last - scalar holds the last element of the right side vector f that
%            is used for least squares solution.
%  nrcount - scalar count used for Arrnorm
%        P - (n x k+1) orthogonal matrix used in the Lanczos bidiagonalization
%        Q - (m x k+1) orthogonal matrix used in the Lanczos bidiagonalization
% reorth012- scalar indicates no, one or two sided reorthogonalization
%       w  - (n x 1) vector used to update soliution x in LSQR
%       x  - (n x 1) vector the LSQR solution
%
% Updated: 1/31/2012
% Authors: James Baglama and Dan Richmond
% Function is called from alsqr
%
%
% First iteration. Set all initial values. 
if k == 0
  beta = norm(Q(:,1)); f(1,1) = beta; Q(:,1) = Q(:,1)/beta;
  if Anumeric 
     P(:,1) = (Q(:,1)'*A)';
  else
     P(:,1) = A(Q(:,1),2);
  end   
  alpha = norm(P(:,1)); P(:,1) = P(:,1)/alpha; B(1,1) = alpha; 
  Arrnorm(1,1)  = alpha*beta;      % Arrnorm(1,1) = ||A'r_0 ||
  Arrnorm(1,2)  =  beta;           % Arrnorm(1,2) = ||r_0||
  Arrnorm(1,3)  = 1;
  nrcount = 1;
end

% This part assumes entry is of the following form
%
%   A * P(:,1:k)   = Q(:,1:k+1)*B(1:k+1,1:k)
%   A'* Q(:,1:k+1) = P(:,1:k)*B((1:k+1,1:k)'+alpha_{k+1}*P(:,k+1)*e_{k+1}'
%
% No checks are performed and the code will produce erroneous results
% if the above equations do not hold on entry. Note k = 0 on first iteration.     

% Set the B(k+1,k+1). On input when k > 0, B(k+1,k+1) holds the alpha value.
alpha = B(k+1,k+1); 

% Compute A*P(:,k+1) - alpha_(k+1)*Q(:,k+1).
if Anumeric
   Q(:,k+2) = A*P(:,k+1)- Q(:,k+1)*alpha;
else
   Q(:,k+2) = A(P(:,k+1),1)- Q(:,k+1)*alpha;
end   
    
% Reorthogonalization of Q vectors.  
if (reorth012 >= 2 || (reorth012 <= 1 && m <= n ) )
    for i=1:k+1
        Q(:,k+2) = Q(:,k+2) - Q(:,i)*(Q(:,i)' * Q(:,k+2));
    end
end
    
% Computes beta_{k+2} of the B(k+2,k+1) matrix
beta = norm(Q(:,k+2)); 
Q(:,k+2) = Q(:,k+2)/beta; 
B(k+2,k+1) = beta;
     
% Compute A'*Q(:,k+2) - beta_{k+2}*P(:,k+1).
if Anumeric
   P(:,k+2) = (Q(:,k+2)'*A)' - P(:,k+1)*beta;
else
   P(:,k+2) = A(Q(:,k+2),2) - P(:,k+1)*beta;
end   
    
% Reorthogonalization of P.
if ((reorth012 >= 2) || (reorth012 <= 1 && m > n ))
   for i=1:k+1
       P(:,k+2) = P(:,k+2) - P(:,i)*(P(:,i)'*P(:,k+2));
   end
end

% B(k+2,k+2) element.
alpha = norm(P(:,k+2));  
B(k+2,k+2) = alpha;
P(:,k+2) = P(:,k+2)/alpha; 
    
% We now have the following relationship
% 
%   A *P(:,1:k+1) = Q(:,1:k+2)*B(1:k+2,1:k+1)
%   A'*Q(:,1:k+2) = P(:,1:k+1)*B((1:k+2,1:k+1)'+alpha_{k+2}*P(:,k+2)*e_{k+2}'
% 
% Right hand side vector is given from the previous iteration and
% is comuputed in the main loop. Notice that 
% f is in the range of Q(:,1:k+1) and Q(:,k+2)'f = 0.
% We need the right hand side to be k+2 x 1 vector.
f(k+2,1) = 0;  f_last=0;
 
% Use MATLAB's internal qr. This gives +- values on the diagonal. 
[Q_B,R_B] = qr(B(1:k+2,1:k+1));
 
% Setup of transition values for LSQR.
rho_bar = alpha*Q_B(k+2,k+2);
theta =   alpha*Q_B(k+2,k+1); 
    
% Solve ||f(k+2,1) - B(k+2,k+1)*y||_2
% Replace B(k+2,k+1) with Q(1:k+2,1:k+2)*R(1:k+2,1:k)
% Replace f(k+2,1) with Q'*f.
f = Q_B'*f; 
phi_bar = f(k+2,1); % Intermediate right hand side value needed for LSQR. 
    
% Compute the solution y of ||f(k+2,1) - B(k+2,k+1)*y||_2.
y = R_B(1:k+1,1:k+1)\f(1:k+1);
  
% Update the solution vector.
x = x + P(:,1:k+1)*y;
    
% Check the stopping criteria.
nrcount = nrcount + 1;
Arrnorm(nrcount,1)  = abs(alpha*beta*y(k+1));%Arrnorm(1,1);      
Arrnorm(nrcount,2)  = abs(phi_bar); %Arrnorm(1,2);           
Arrnorm(nrcount,3)  = Arrnorm(nrcount-1,3) + 2;
if Arrnorm(nrcount,1) <= ltol *Arrnorm(1,1)
   flag(1) = 1; Arrnorm = Arrnorm(1:nrcount,1:3); return;
end
    
% Update w vector
w = P(:,k+2) - P(:,1:k+1)*(y/f(k+1,1)*theta);
    
% Setup that last vectors to hold the current iteration value
% for running the LSQR iteration.
P(:,m_b+1) = P(:,k+2); Q(:,m_b+1) = Q(:,k+2);
    
for i = k+2:maxitl
    
   %--------------------------------------%
   % Sart Lanczos Bidiagonalization steps %
   %--------------------------------------%
     
   % Set the next alpha value to diagonal value of B(i,i).
   if i <= m_b % Building B matrix.
        B(i,i) = alpha; 
   end
     
   % Compute A*v - alpha*u.
   if Anumeric
      Q(:,m_b+1) = A*P(:,m_b+1) - alpha*Q(:,m_b+1); 
   else
      Q(:,m_b+1) = A(P(:,m_b+1),1) - alpha*Q(:,m_b+1);
   end   
      
   % Reorthogonalization of Q vectors.
   if i <= m_b
     if (reorth012 >= 2 || (reorth012 <= 1 && m <= n ) )
          for j=1:i
             Q(:,m_b+1) = Q(:,m_b+1) - (Q(:,j)'*Q(:,m_b+1))*Q(:,j);
          end
   elseif(reorth012 >= 2 || (reorth012 == 1 && m <= n ) )
       for j=1:min(i,m_b)
           Q(:,m_b+1) = Q(:,m_b+1) - (Q(:,j)'*Q(:,m_b+1))*Q(:,j); 
       end
     end
   end
    
   % Computes beta_{i} of the B matrix
   beta = norm(Q(:,m_b+1)); Q(:,m_b+1) = Q(:,m_b+1)/beta;
     
   % Set the beta value to B(i+1,i).
   if i <= m_b, B(i+1,i) = beta; end % Building B matrix.
     
   % Compute A'*Q - beta*P.
   if Anumeric
      P(:,m_b+1) = (Q(:,m_b+1)'*A)' - beta*P(:,m_b+1);
   else
      P(:,m_b+1) = A(Q(:,m_b+1),2) - beta*P(:,m_b+1);
   end   
    
   % Reorthogonalization of P vectors.
   if i <= m_b
     if (reorth012 >= 2 || (reorth012 <= 1 && m > n ) )
         for j=1:i
             P(:,m_b+1) = P(:,m_b+1) - (P(:,j)'*P(:,m_b+1))*P(:,j);
         end   
     elseif(reorth012 >= 2 || (reorth012 == 1 && m > n ) )
          for j=1:min(i,m_b)
              P(:,m_b+1) = P(:,m_b+1) - (P(:,j)'*P(:,m_b+1))*P(:,j);
          end
     end
   end
 
   % Computes alpha_{i} of the B matrix    
   alpha = norm(P(:,m_b+1)); P(:,m_b+1) = P(:,m_b+1)/alpha;
    
   %-------------------------------------%
   % End Lanczos Bidiagonalization steps %
   %-------------------------------------%

   % Store only Q(:,1:m_b) and P(:,1:m_b) 
   % and the matrix B(1:m_b+1,1:m_b) 
   % Storage of the Q and P vectors is overwritten.
   % and reorthogonalization occcurs with all k harmonic Ritz
   % vectors and the last generated m_b vectors.
   if i < m_b
        Q(:,i+1) = Q(:,m_b+1); P(:,i+1) = P(:,m_b+1);
   elseif m_b ~= maxitl
       Q(:,mod(i,m_b)+1) = Q(:,m_b+1); 
       P(:,mod(i,m_b)+1) = P(:,m_b+1);   
   end
   
    %-----------------------------%
    % Start LSQR solution process %
    %-----------------------------%
    
    % Givens values. 
    rho = pythag(rho_bar,beta); c1 = rho_bar/rho; s1 = beta/rho; 
    
    % Translate B from lower bidiagonal to upper bidiagonal matrix R.
    % Super diagonal element in R is theta_{i} and is computed
    % from applying Givens to B. The diagonal value of R is rho_{i}
    % and is given above as two norm of [rho_bar,beta].
    theta = s1*alpha; 
    
    % Updating the last diagonal element of R.
    rho_bar = -c1*alpha;
    
    % Updating elements in f (two elements up).
    phi = c1*phi_bar; phi_bar = s1*phi_bar;
    f_last = -c1*phi_bar;

    % Update the solution.
    x = x + (phi/rho)*w; w = P(:,m_b+1) - (theta/rho)*w;

    %---------------------------%
    % End LSQR solution process %
    %---------------------------% 
    
    % Calculate the residual norms ||B-A*X(J)||
    % and ||A'(*B- A*X(J))|| and compute
    % matrix-vector products used to compute X(J).
    % ||A'*(B - A*X)|| <= TOLLSQR*||A||*||B - A*X||
    % Check the stopping criteria.
    nrcount = nrcount + 1; flag(4) = i;
    Arrnorm(nrcount,1)  = abs(phi_bar*rho_bar);%Arrnorm(1,1);      
    Arrnorm(nrcount,2)  = abs(phi_bar); %Arrnorm(1,2);           
    Arrnorm(nrcount,3)  = Arrnorm(nrcount-1,3) + 2;
    if Arrnorm(nrcount,1) <= ltol*Arrnorm(1,1)
       flag(1) = 1;  
       Arrnorm = Arrnorm(1:nrcount,1:3); 
       return;
    end
end 