function [v,theta,acc,estconv,trueconv,allconv,backerr] = ...
             IDRRitzPurified(A,R,uSH,SR,PR,omega,resvec,lambda,s,maxj)

%IDRRITZPURIFIED is a howto: estimate accuracy without long vectors.
%
%   [v,theta,acc,estconv,trueconv,allconv,backerr] = ...
%         IDRRitzPurified(A,R,uSH,SR,PR,omega,resvec,lambda,s,maxj)
%
%   computes with given
%
%    A        - the square system matrix,
%    R        - the matrix of all IDRORes residuals,
%    (uSH,SR) - extended Sonneveld pencil,
%    (uSH,PR) - extended purified pencil,
%    omega    - vector of all chosen omega,
%    resvec   - vector of norms of all IDRORes residuals,
%    lambda   - eigenvalue of A of interest,
%    s        - parameter `s' of IDR(s)ORes,  and
%    maxj     - number of sweeps of IDRORes
%
%  the values
%
%    v - normalized eigenvector of the purified pencil,
%    theta - corresponding Ritz value,
%    acc - final convergence estimator,
%    estconv - vector of convergence estimators of all steps,
%    trueconv - vector of smallest distances to lambda of all steps,
%    allconv - matrix of all distances to lambda of all steps,  and
%    backerr - vector of backward errors of best Ritz pairs of all
%      steps.
%
%   See also: IDRORes, IDRORes_academic, IDREig, IDREig_driver.

%   Coypright 2009-2010 Jens-Peter M. Zemke.
%   Version 1.01, Date: 2010/04/13 14:57:36 CEST.

n = (s+1)*maxj;
omega_theta = ones(n,1);
allconv = zeros(n);
estconv = zeros(n,1);
trueconv = zeros(n,1);
backerr = zeros(n,1);

for i = 1:n

  %%% compute the eigenpairs of the purified pencil
  [S,Theta] = eig(full(uSH(1:i,1:i)),full(PR(1:i,1:i)));
  dtheta = diag(Theta).';

  %%% sort according to the distance to the given lambda
  [allconv(i,1:i),index] = sort(abs(lambda-dtheta));
  theta = dtheta(index(1));
  trueconv(i) = min(allconv(i,1:i));

  %%% compute the scaling by polynomial terms
  for k = s+1:i
    omega_theta(k) = prod(1-(omega(1:floor(k/(s+1))).*theta));
  end

  %%% compute the approximate eigenvector of the Sonneveld pencil
  v = S(:,index(1))./omega_theta(1:i);

  %%% still cheap: refine using Sonneveld upper triangular matrix
  w = SR(1:i,1:i)*v;

  %%% costly: prolongation using the matrix of residual vectors
  w = R(:,1:i)*w;

  %%% costly: computation of the backward error
  w = w/norm(w);
  backerr(i) = norm(A*w-theta*w);

  %%% normalize the short vector at hand
  v = v/norm(v);

  %%% estimate cheaply the rate of convergence
  estconv(i) = resvec(i+1)*abs(uSH(i+1,i)*v(i));

end

acc = estconv(n);
