function [RW,pnRW,W,pnW] = compute_W(A,L,r,omega,resvec,s,maxj)

%COMPUTE_W computes the basis vectors in W for IDRORes.
%
%   [RW,pnRW,W,pnW] = compute_W(A,L,r,omega,resvec,s,maxj)
%   computes the matrix W consisting column-wise of the basis
%   vectors of the BiORes(s,1)-process underlying IDR(s)ORes.
%   This is done using the Lanczos-BiORes(s,1) matrix computed
%   from IDREig with the information obtained and stored by
%   IDRORes_academic.
%
%   On input,
%
%      A      - denotes the quadratic system matrix,
%      L      - denotes the Lanczos-BiORes(s,1) matrix
%               computed by IDREig,
%      r      - is the starting residual used for IDRORes,
%      omega  - is the vector of all omega chosen in IDRORes,
%      resvec - contains the norms of the updated residuals,
%      s      - denotes the parameter `s' of IDR(s)ORes, and
%      maxj   - denotes the number of sweeps carried out.
%
%   On output,
%
%      RW   - contains the recomputed residuals r,
%      pnRW - is a prolonged vector of norms of these recomputed
%             residuals,
%      W    - contains the computed Lanczos-BiORes(s,1) residual
%             vectors, and
%      pnW  - is a prolonged vector of norms of these computed
%             Lanczos-BiORes(s,1) residuals.
%
%   See also: IDRORes_academic, IDREig, IDREig_driver.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                         %%
%%   use the OrthoRes-Lanczos-matrix for the computation   %%
%%              of the residual polynomials                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m = maxj*s;
n = size(A,1);
maxit = (s+1)*maxj;

W = zeros(n,m+1);

W(:,1) = r;
nW = eye(m,1)*norm(W(:,1));

for j = 1:m-1
  W(:,j+1) = (A*W(:,j)-W(:,1:j)*L(1:j,j))/L(j+1,j);
  nW(j+1) = norm(W(:,j+1));
end

indvec = ((1:maxit) ~= (s+1)*floor((1:maxit)./(s+1))) ~= 0;

%% prolong nW

pnW = zeros(maxit,1);
pnW(indvec) = nW;

%%% compute R by Omega_j(A)*W by a Horner scheme

RW = W;
nRW = nW;

for j = 1:maxj-1
  RW(:,j*s+1:m) = RW(:,j*s+1:m)-omega(j)*A*RW(:,j*s+1:m);
  for k = 1:s
    nRW(j*s+k) = norm(RW(:,j*s+k));
  end
end

pnRW = zeros(maxit,1);
pnRW(indvec) = nRW;
