function [U,H,flag] = arnoldi_W(A,ns,x)
% Realiza ns iteraciones de Arnoldi con reortogonalizavión en A partiendo de x
% Si flag == 0, las columans de U are ortonormales, H(ns+1,ns) Hessenberg matrix y
%    A*U(:,1:ns) = U(:,1:ns+1)*H.
% Si flag == j > 0, se ha obtenido subespacio invariante en j iter. y
%    A*U(:,1:j) = U(:,1:j)*H(1:j,1:j).
%
 flag = 0;
 H = zeros(ns+1,ns); n      = size(A,1);
 U = zeros(n,ns+1);  U(:,1) = x/norm(x); delta = zeros(ns+1,1);
 for jj = 1:ns
    U(:,jj+1)   = A*U(:,jj);                         % sólo se usa A aquí
    H(1:jj,jj)  = U(:,1:jj)'*U(:,jj+1);
    U(:,jj+1)   = U(:,jj+1) - U(:,1:jj)*H(1:jj,jj);  % ortogonaliza
    delta(1:jj) = U(:,1:jj)'*U(:,jj+1);
    U(:,jj+1)   = U(:,jj+1) - U(:,1:jj)*delta(1:jj); % reortogonaliza
    H(1:jj,jj)  = H(1:jj,jj) + delta(1:jj);
    H(jj+1,jj)  = norm(U(:,jj+1));
    if H(jj+1,jj) == 0, flag = jj; return, end
    U(:,jj+1)   = U(:,jj+1)/H(jj+1,jj);
 end
end