function [uSH,SR,uSM,PR,uDH,DR,uLM,D] = IDREig(uH,UT,C,omega,s,maxj);

%IDREIG constructs matrices for eigenvalue computations for IDRORes.
%
%   [uSH,SR] = IDREig(uH,UT,C,omega,s,maxj) computes the sparse
%   extended Sonneveld pencil from the quantities returned by a
%   run of the IDRORes implementation of the prototype IDR(s) of
%   [Sonneveld/van Gijzen, SISC 31, 2008].
%
%   In the call
%
%   [uSH,SR,uSM,PR,uDH,DR,uLM] = IDREig(uH,UT,C,omega,s,maxj)
%
%     - uH is the leading extended Hessenberg matrix defined by
%       the starting procedure,
%     - UT is the upper triangular matrix defined by the starting
%       procedure,
%     - C contains the coefficients of the oblique projections
%       onto the differences of the last residuals to make the
%       resulting vector orthogonal to the shadow vectors in P,
%     - omega contains the vector of all chosen omega,
%     - s is a non-negative integer which defines the size of
%       the shadow space,  and
%     - maxj gives the number of sweeps.
%
%  From the output,
%
%     (uSH,SR) is the sparse extended Sonneveld pencil,
%      uSM     is the sparse extended Sonneveld matrix,
%     (uSH,PR) is the sparse extended purified pencil,
%     (uDH,DR) is the sparse extended deflated pencil,  and
%      uLM     is the sparse extended Lanczos-BiORes(s,1) matrix.
%
%   See also: IDRORes, IDRORes_academic, IDREig_driver, IDRRitzPurified.

%   Coypright 2009-2010 Jens-Peter M. Zemke.
%   Version 1.00, Date: 2010/04/13 14:33:02 CEST.

j = maxj-1;
n = (s+1)*j+1;

o = zeros(1,n);
e = ones(1,n);

band_Yo = diff([o;C;e;o]);
band_Y  = diff([o;C;e]);

diag_omega = diag(kron(diag(sparse([1,omega(1:j)])),speye(s+1)));
D = spdiags([diag_omega(2:end);omega(end)],0,n+s,n+s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                   %%
%%        compute the extended Sonneveld pencil (uSH,SR)             %%
%%                                                                   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uSH = [[sparse(uH);sparse(n,s)],...
        spdiags(fliplr(band_Yo.'),[-(s+1):0],n+s+1,n)];
SR  = [[sparse(UT);sparse(n,s)],...
        spdiags(fliplr(band_Y.'), [-s:0],n+s,n)]*D;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                   %%
%%        compute the extended Sonneveld matrix uSM                  %%
%%                                                                   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uSM = uSH/SR;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                   %%
%%       compute the matrix PR of the purified pencil                %%
%%                                                                   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

band_U = band_Y;
band_U(repmat(fliplr(triu(logical(ones(s+1)))),1,j)) = 0;
band_U(:,n) = 0;
PR = [[sparse(UT);sparse(n,s)],...
       spdiags(fliplr(band_U.'),[-s:0],n+s,n)]*D;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                   %%
%%       compute the deflated purified pencil (uDH,DR)               %%
%%                                                                   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G = speye(n+s);
indicesG = (1:maxj).*(s+1);
for refine = 1:2
  G(indicesG,:) = -uSH(indicesG,indicesG)\uSH(indicesG,:);
  G(indicesG,indicesG) = speye(maxj);
  uSHG = uSH*G;
end
indicesD = (1:n+s) ~= (s+1)*floor((1:n+s)./(s+1));

uDH = uSHG([indicesD logical(1)],indicesD);
DR = PR(indicesD,indicesD);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                   %%
%%       compute the BiORes(s,1) matrix uLM                          %%
%%                                                                   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uLM = uDH/DR;
