function [iter,V,D] = albeigs(atv,n,blsz,K,tol,varargin)
%
% find a few eigenvalues and eigenvectors for the standard eigenvalue problem
% A*x = lambda*x using the implicit restarted Arnoldi method (IRAM).
% based on a Matlab script by James Baglama, University of Rhode Island, 2005.
% [iter,V,D] = albeigs(atv,n,blsz,K,tol,varargin)
% atv    character name of a matrix-vector product M-file returning Ab where b
%        is input. The format for atv is "function x=atv(b,n,blsz,iter,p1,p2,...)"
%        where p1 and p2 are optional parameters
% n      size of matrix A
% blsz   block size of the Arnoldi Hessenberg matrix (blsz=3 recommended)
% K      number of desired eigenvalues
% tol    tolerance used for convergence (tol=1.0d-6 recommended)
% varargin optional parameters (p1,p2,...) for atv
% output parameters:
% iter   number of Arnoldi iterations.
% V      (N x blsz) eigenvector matrix
% D      (N x blsz) diagonal eigenvalue matrix
% (c) 2020 Alain Hebert, Polytechnique Montreal

  maxit = 100; nbls = 10;

  % Resize Krylov subspace if blsz*nbls (i.e. number of Arnoldi vectors) is
  % larger than n (i.e. the size of the matrix A).
  if blsz*nbls >= n
    nbls = floor(n/blsz);
    warning(['Changing nbls to ',num2str(nbls)]);
  end

  % Increase the number of desired values to help increase convergence.
  % Set K+adjust to be next multiple of blsz > K.
  vadjust = 1:blsz; I = find(mod(K + vadjust-1,blsz) == 0); adjust = I(1);
  K_org = K; K = K + adjust; % K_org is the orginal value of K.

  % Check for input errors in the structure array.
  if K <= 0,    error('K must be a positive value.'),        end
  if K > n,     error('K is too large. Use eig(full(A))!'),  end
  if blsz <= 0, error('blsz must be a positive value.'),     end
  if tol < 0,   error('TOL must be non-negative.'),          end

  % Automatically adjust Krylov subspace to accommodate larger values of K.
  if blsz*(nbls-1) - blsz - K - 1 < 0
    nbls = ceil((K+1)/blsz+2.1);
    warning(['Changing nbls to ',num2str(nbls)]);
  end
  if blsz*nbls >= n
    nbls = floor(n/blsz-0.1);
    warning(['Changing nbls to ',num2str(nbls)]);
  end
  if blsz*(nbls-1) - blsz - K - 1 < 0,  error('K is too large.'); end
  V = [tril(ones(n,blsz)), zeros(n,nbls*blsz)];
  if tol < eps, tol = eps; end % Set tolerance to machine precision if tol < eps.
  iter = 0;     % Main loop iteration count.
  m = blsz;     % Current number of columns in the matrix V.
  T = [];       % T matrix in the WY representation of the Householder products.

  while 1
    iter = iter+1;
    if iter > maxit
      error('maximum number of IRAM iterations reached')
    end
    % Compute the block Householder Arnoldi decomposition.
    T_blk =zeros(blsz,blsz);
    if blsz*(nbls+1) > n, nbls = floor(n/blsz)-1; end
    % Main iteration loop for the Augmented Block Householder Arnoldi decomposition.
    while (m <= blsz*(nbls+1))
      [V(m-blsz+1:n,m-blsz+1:m),tau] = alst2c(real(V(m-blsz+1:n,m-blsz+1:m)));
      if m > blsz
        H(1:m-blsz,m-2*blsz+1:m-blsz) = V(1:m-blsz,m-blsz+1:m);
        H(m-blsz+1:m,m-2*blsz+1:m-blsz) = triu(V(m-blsz+1:m,m-blsz+1:m));
      end
      m2 = 2*m;
      if m2 > n, m2=n; end
      V(1:m2,m-blsz+1:m) = tril(V(1:m2,m-blsz+1:m),-(m-blsz+1));
      for i=m-blsz+1:m, V(i,i)=1; end
      T_blk(1,1) = -2/(V(:,m-blsz+1)'*V(:,m-blsz+1));
      for i=2:blsz
        T_blk(1:i,i) = (V(:,m-blsz+i)'*V(:,m-blsz+1:m-blsz+i))';
        T_blk(i,i) = -2/T_blk(i,i);
        T_blk(1:i-1,i) = T_blk(i,i)*T_blk(1:i-1,1:i-1)*T_blk(1:i-1,i);
      end
      T = [ [T; zeros(blsz,m-blsz)] ...
            [(T*(V(:,m-blsz+1:m)'*V(:,1:m-blsz))')*T_blk; T_blk] ];
      if m <= blsz*nbls % Resize and reactualize V
        V(:,m+1:m+blsz) = V(:,1:m)*(T*V(m-blsz+1:m,1:m)');
        V(m-blsz+1:m,m+1:m+blsz) = eye(blsz,blsz)+V(m-blsz+1:m,m+1:m+blsz);
        V(:,m+1:m+blsz) = feval(atv,V(:,m+1:m+blsz),n,blsz,iter,varargin{:});
        V(:,m+1:m+blsz) = V(:,m+1:m+blsz) + V(:,1:m)*((V(:,m+1:m+blsz)'*V(:,1:m))*T)';
      end
      m = m + blsz;
    end

    % Determine the size of the block Hessenberg matrix H. Possible truncation
    % may occur if an invariant subspace has been found.
    H=real(H);
    [Hsz_n,Hsz_m] = size(H);

    % Compute the eigenvalue decomposition of the block Hessenberg H(1:Hsz_m,:).
    [jter,H_eigv,H_eig] = alhqr(H(1:Hsz_m,:),200);
    H_eig = diag(H_eig);

    % Check the accuracy of the computation of the eigenvalues of the
    % Hessenberg matrix. This is used to monitor balancing.
    conv = 'F';  % Boolean to determine if all desired eigenpairs have converged.

    % Sort the eigenvalue and eigenvector arrays.
    [Y,I_eig] = sort(abs(H_eig)); I_eig = I_eig(Hsz_m:-1:1);
    H_eig = H_eig(I_eig); H_eigv = H_eigv(:,I_eig);

    % Compute the residuals for the K_org Ritz values.
    residuals = (sum(abs(H(Hsz_n-blsz+1:Hsz_n,Hsz_m-blsz+1:Hsz_m)*...
                         H_eigv(Hsz_m-blsz+1:Hsz_m,1:K_org)).^2, 1).^(0.5));

    % Check for convergence.
    Len_res = length(find(residuals(1:K_org)' < tol*abs(H_eig(1:K_org))));
    if Len_res >= K_org, conv = 'T'; end % Set convergence to true.

    % Adjust K to include more vectors as the number of vectors converge.
    Len_res = length(find(residuals(1:K_org)' < eps*abs(H_eig(1:K_org))));
    K =  K_org + adjust+ Len_res; if K > Hsz_m - 2*blsz-1, K = Hsz_m - 2*blsz-1; end

    % Determine if K splits a conjugate pair. If so replace K with K + 1.
    if ~isreal(H_eig(K))
      eig_conj = 1;
      if K < Hsz_m
        if abs(imag(H_eig(K)) + imag(H_eig(K+1))) < sqrt(eps)
          K = K + 1; eig_conj = 0;
        end
      end
      if K > 1 && eig_conj
        if abs(imag(H_eig(K)) + imag(H_eig(K-1))) < sqrt(eps)
          eig_conj = 0;
        end
      end
      if eig_conj, error('%d th conjugate pair split.',K); end
    end

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

    % Compute the QR factorization of H_eigv(:,1:K).
    if ~isreal(H_eig)
      % Convert the complex eigenvectors of the eigenvalue decomposition of H
      % to real vectors and convert the complex diagonal matrix to block diagonal.
      iset = find(imag(H_eig)'); index = iset(1:2:end);
      if isempty(index)
        Q = H_eigv(:,1:K);
      else
        if max(index)==size(H_eig,1) || any(conj(H_eig(index))~=H_eig(index+1))
          error('invalid diagonal');
        end
        t = eye(size(H_eig,1));
        twobytwo = [1 1;sqrt(-1) -sqrt(-1)];
        for i=index
          t(i:i+1,i:i+1) = twobytwo;
        end
        Q=H_eigv/t;
      end
    else
      Q = H_eigv(:,1:K);
    end

    % Economy size QR factorization
    [X, tau] = alst2c(real(Q(:,1:K)));
    R = triu(X(1:K,1:K)); Q = eye(Hsz_m,K);
    for j = K:-1:1
      w = [1; X(j+1:Hsz_m,j)];
      Q(j:Hsz_m,:) = Q(j:Hsz_m,:)+tau(j)*w*(w'*Q(j:Hsz_m,:));
    end

    % The Schur matrix for H.
    T_Schur = triu(Q'*H(1:Hsz_m,:)*Q,-1);

    % Compute the starting vectors and the residual vectors from the Householder
    % WY form. The starting vectors will be the first K Schur vectors and the
    % residual vectors are stored as the last blsz vectors in the Householder WY form.
    Tsz_m = size(T,1);
    V(:,Hsz_n-blsz+1:Hsz_n)= V(:,1:Tsz_m)*(T*V(Tsz_m-blsz+1:Tsz_m,1:Tsz_m)');
    V(Tsz_m-blsz+1:Tsz_m,Hsz_n-blsz+1:Hsz_n) = eye(blsz,blsz)+ ...
        V(Tsz_m-blsz+1:Tsz_m,Hsz_n-blsz+1:Hsz_n);
    V(:,1:K) = V(:,1:Hsz_m)*(T(1:Hsz_m,1:Hsz_m)*(V(1:Hsz_m,1:Hsz_m)'*Q(:,1:K)));
    V(1:Hsz_m,1:K) = Q(:,1:K) + V(1:Hsz_m,1:K);

    % Set the size of the large matrix V and move the residual vectors.
    m = K + 2*blsz; V(:,K+1:K+blsz) = V(:,Hsz_n-blsz+1:Hsz_n);

    % Set the new starting vector(s) to be the desired vectors V(:,1:K) with the
    % residual vectors V(:,Hsz_n-blsz+1:Hsz_n). Place all vectors in the compact
    % WY form of the Householder product. Compute the next set of vectors by
    % computing A*V(:,Hsz_n-blsz+1:Hsz_n) and store this in V(:,Hsz_n+1:Hsz_n+blsz).
    m2=m-blsz;
    V_sign = zeros(m2,1);
    D = zeros(m2,m2); T = V(1:m2,1:m2);
    V(:,m2+1:m) = feval(atv,V(:,m2-blsz+1:m2),n,blsz,iter,varargin{:});
    for i =1:m2
      V_sign(i) = sign(V(i,i));
      if V_sign(i) == 0, V_sign(i)=1; end
      Vdot = 1 + V_sign(i)*V(i,i);      % Dot product of Householder vectors.
      V(i,i) = V(i,i) + V_sign(i);      % Reflection to the ith axis.
      D(i,i) = 1/V(i,i);                % Used for scaling. Note: V(i,i) >= 1.
      V(1:m2,i+1:m2) = V(1:m2,i+1:m2) - V(1:m2,i)*((V_sign(i)/Vdot)*V(i,i+1:m2));
    end
    V(1:m2,1:m2) = V(1:m2,1:m2)*D;
    R = -(diag(V_sign));
    V(:,m2+1:m) = V(:,m2+1:m)*R(m2-blsz+1:m2,m2-blsz+1:m2);
    V(1:m2,1:m2) = tril(V(1:m2,1:m2),0);
    T = T*R - eye(m2); T = T/V(1:m2,1:m2)';  T = triu(V(1:m2,1:m2)\T);
    V(m2+1:n,1:m2) = V(m2+1:n,1:m2)*(R/(T*V(1:m2,1:m2)'));
    V(:,m2+1:m) = V(:,m2+1:m) + V(:,1:m2)*((V(:,m2+1:m)'*V(:,1:m2))*T)';

    % Compute the first K columns and K+blsz rows of the matrix H, used in augmenting.
    H_R = H(Hsz_n-blsz+1:Hsz_n,Hsz_m-blsz+1:Hsz_m);
    H=zeros(K+blsz,K); % Resize H
    H(1:K,1:K) = R(1:K,1:K)*T_Schur(1:K,1:K)*R(1:K,1:K);
    H(K+1:K+blsz,1:K) = R(K+1:K+blsz,K+1:K+blsz)*H_R* ...
                        Q(Hsz_m-blsz+1:Hsz_m,1:K)*R(1:K,1:K);
  end

  % Truncated eigenvalue and eigenvector arrays to include only desired eigenpairs.
  H_eig = H_eig(1:K_org); H_eigv = H_eigv(:,1:K_org);
  Tsz_m = size(T,1); [H_sz1,H_sz2] = size(H_eigv);
  V = V(:,1:Tsz_m)*(T*(V(1:H_sz1,1:Tsz_m)'*H_eigv));
  V(1:H_sz1,1:H_sz2) = H_eigv + V(1:H_sz1,1:H_sz2);
  D = diag(H_eig);