function [H,V] = arn_pl(A, E, k,w,opts)
%
%  Arnoldi method w.r.t. inv(E)*A
%
%  Calling sequence:
%
%    [H,V] = arn_p(A, E, k,w,opts)

%  Input:
%
%    A,E         the matrices A;
%    k         number of Arnoldi steps (usually k << n, where
%              n is the dim. of A);
%    w         initial n-vector 
%               (optional - chosen at random, if omitted).
%    opts      optional structure contains settings for system solve 
%               (if not present, direct solves are used)    
%  Output:
%
%    H         matrix H, (k+1)-x-k matrix, upper Hessenberg;
%    V         matrix V, n-x-(k+1) matrix, (nearly) orthogonal columns,
% 
% so that V(:,1) in span{w},  V'*V = eye(k+1), inv(E)*A*V(:,1:k) = V*H.
%   This might hold only approximately, if the linear systems are solved
%   approximately

% Input data is not completely checked, nor does this code check for
% breakdown.
   

na = nargin;
n = size(A,1);                 
if k >= n-1, error('k must be smaller than the order of A!'); end
if na<4, w = randn(n,1); end
if na<5, opts=[]; end

V = zeros(n,k+1);
H = zeros(k+1,k);

V(:,1) = (1.0/norm(w))*w;

beta = 0;

for j = 1:k
  if j > 1
    H(j,j-1) = beta;
    V(:,j) = (1.0/beta)*w;
  end
  
  w = MV(A,V(:,j));
  if ~isempty(E) % direct solve (default)
      if isempty(opts)
          w=E\w;
      else % solve as specified by opts 
          w=solve_LS(@(x) MV(E,x),w, opts.inner_tol,opts.maxit_inner,opts.linsolver,[],opts.M1,opts.M2,opts.Mtype);
      end
  end
  for it=1:2 % double orthogonalization (is enough)
      for i = 1:j
          coef = V(:,i)'*w;
          H(i,j) = H(i,j)+coef;
          w = w-coef*V(:,i);
      end
  end
  beta = norm(w);
  H(j+1,j) = beta;
end  
V(:,k+1) = (1.0/beta)*w;