function [x,resvec,uH,UT,C,omega] = IDRORes(A,b,x0,s,PH,maxj,start)

%IDRORes approximate solutions to Ax=b using original IDR.
%
%   x = IDRORes(A,b) computes an approximate solution x to the
%   linear system Ax = b using the IDR(s) prototype algorithm
%   by [Sonneveld/van Gijzen, SISC 31, 2008].
%
%   In the call
%   [x,resvec,uH,UT,C,omega] =
%                          IDRORes(A,b,x0,s,PH,maxj,start)
%   of IDR(s)ORes,
%  
%     - x0    is the starting guess,
%     - s     is the number of left Lanczos vectors,
%     - PH    is the left Petrov-Galerkin matrix,
%     - maxj  is the number of inner sweeps
%                (comprising of s left and 1 right steps),
%     - start is one of:
%          1 == OrthoRes(1) = ORes,
%          2 == full OrthoRes,
%
%   and
%
%     - x       is the approximate solution,
%     - resvec  is the vector of the norms of the
%                  computed residuals,
%     - uH      is a sparse unreduced extended upper
%                  Hessenberg matrix,
%     - UT      is a sparse upper triangular matrix,
%     - C       is a (s x n)-matrix that contains in the
%                  columns the coefficients obtained in
%                  the orthogonalization against P,
%     - omega   contains the local minimizers.
%
%   The extended pencil (uH,UT) is defined by the starting
%   procedure 1 or 2. This pencil, together with the IDR
%   quantities contained in the matrix C and the vector
%   omega, defines the underlying Hessenberg decomposition
%   for the residuals that can be used to compute eigenvalue
%   approximations.
%
%   There is also an academic implementation with two drawbacks:
%     o The computation of the true residual norm is based on
%       a second matrix-vector multiplication in every step,
%     o All the residuals are stored in the matrix R.
%   The academic variant is named IDRORes_academic.
%
%   See also: IDRORes_academic, IDREig, IDREig_driver, IDRRitzPurified.

%   Coypright 2009-2010 Jens-Peter M. Zemke.
%   Version 1.03, Date: 2010/03/10 16:31:43 CEST.

%% TODO: implement full GMRes

if nargin < 1
  error('No system matrix was given.');
end

[n,m] = size(A);

if nargin < 2
  error('No right-hand side was given.');
end

if nargin < 3
  x0 = randn(n,1);
end

if nargin < 4
  s = min(n,4);
end

if nargin < 5
  PH = orth(randn(n,s)).';
end

if nargin < 6
  maxj = min(round((2*(n+1))/(s+1)),round(300/(s+1)));
end

if nargin < 7
  start = 1;
end

resvec = zeros(1,maxj*(s+1)+1);

x = x0;
r = b-A*x0;

resvec(1) = norm(r);

UT = sparse(s);
uH = sparse(s+1,s);

  nablaR = zeros(n,s);
  nablaX = zeros(n,s);
PHnablaR = zeros(s,s);

% ORes, i.e., OrthoRes(1)
if start == 1

  for k = 1:s

    Ar = A*r;
    tomega = dot(Ar,r)/dot(Ar,Ar);
    UT(k,k) = tomega;

    nablaX(:,k) = tomega*r;
    x = x + nablaX(:,k);

    nablaR(:,k) = -tomega*Ar;
    r = r - tomega*Ar;

    PHnablaR(:,k) = PH*nablaR(:,k);

    resvec(k+1) = norm(r);

  end
  uH = spdiags(ones(s+1,1)*[-1 1],[-1 0],s+1,s);

% full OrthoRes, i.e., OrthoRes(s+1)
elseif start == 2

  R = zeros(n,s+1);
  X = zeros(n,s+1);
  rHr = zeros(s+1,1);

  R(:,1) = r;
  X(:,1) = x;
  rHr(1) = dot(r,r);

  for k = 1:s

    Ar = A*R(:,k);

    uH(1:k,k) = (R(:,1:k)'*Ar)./rHr(1:k);
    uH(k+1,k) = -sum(uH(1:k,k));

    R(:,k+1) =  (Ar-R(:,1:k)*uH(1:k,k))/uH(k+1,k);
    X(:,k+1) = -(R(:,k)+X(:,1:k)*uH(1:k,k))/uH(k+1,k);

    resvec(k+1) = norm(R(:,k+1));
    rHr(k+1) = resvec(k+1)^2;

  end

  UT = speye(s);
  r = R(:,s+1);
  x = X(:,s+1);

  nablaR = R(:,2:s+1)-R(:,1:s);
  nablaX = X(:,2:s+1)-X(:,1:s);
  PHnablaR = PH*nablaR;

% other methods, e.g., full explicit GMRes, i.e., GMRes(s+1)
else

    disp('not implemented yet')

end

PHr = PH*r;
normr = norm(r);

sp1 = s+1;
step = sp1+1;
j = 1;
colnum = 1;

indexR = (1:s)';
indexC = (1:s)';

C = zeros(s,(maxj-1)*(s+1));
omega = zeros(1,maxj);

for j = 1:maxj-1

  c = (PHnablaR)\(PHr);
  C(:,colnum) = c(indexC);

  oldest = indexR == 1;
  indexR = circshift(indexR, 1);
  indexC = circshift(indexC,-1);

  update = -nablaR*c;
  v = r + update;
  Av = A*v;
  omega(j) = dot(Av,v)/dot(Av,Av);

  nablaX(:,oldest) = -nablaX*c + omega(j)*v;
  x = x + nablaX(:,oldest);

  nablaR(:,oldest) = update - omega(j)*Av;
  r = r + nablaR(:,oldest);

  PHnablaR(:,oldest) = PH*nablaR(:,oldest);
  PHr = PHr + PHnablaR(:,oldest);

  normr = norm(r);
  resvec(step) = normr;

  for k = 1:s

    c = (PHnablaR)\(PHr);
    C(:,colnum+k) = c(indexC);

    oldest = indexR == 1;
    indexR = circshift(indexR, 1);
    indexC = circshift(indexC,-1);

    update = -nablaR*c;
    v = r + update;

    nablaX(:,oldest) = -nablaX*c + omega(j)*v;
    x = x + nablaX(:,oldest);

    nablaR(:,oldest) = -A*nablaX(:,oldest);
    r = r + nablaR(:,oldest);

    PHnablaR(:,oldest) = PH*nablaR(:,oldest);
    PHr = PHr + PHnablaR(:,oldest);

    normr = norm(r);
    resvec(step+k) = normr;

  end

  colnum = colnum + sp1;
  step = step + sp1;

end

c = (PHnablaR)\(PHr);
C(:,colnum) = c(indexC);

oldest = indexR == 1;
indexC = circshift(indexC,-1);

update = -nablaR*c;
v = r + update;
Av = A*v;
omega(j+1) = dot(Av,v)/dot(Av,Av);

nablaX(:,oldest) = -nablaX*c + omega(j+1)*v;
x = x + nablaX(:,oldest);

nablaR(:,oldest) = update - omega(j+1)*Av;
r = r + nablaR(:,oldest);

normr = norm(r);
resvec(step) = normr;
