clear

format long

A=[-2 1 2 4; 1 0 -1 2;-3 2 -1 -1;0 1 -2 1];
b=[1 1 1 1]';
[m,n]=size(A);

x0=[1 0 -1 0]';
r0=b-A*x0;
rho=norm(r0);

h=zeros(n+1,n);
q=zeros(n,n);
q(:,1)=r0/rho;

e1=zeros(n,1);e1(1)=1;

for k=1:n
    for p=1:k
        h(p,k)=q(:,p)'*A*q(:,k);
    end
    t=A*q(:,k);
    for p=1:k
        t=t-h(p,k)*q(:,p);
    end
    h(k+1,k)=norm(t);
    if k<n
        q(:,k+1)=t/h(k+1,k);
    end
    [Q1,R1]=qr(h(1:1+k,1:k),0);
    y=rho*(R1\(Q1(1,:)'));
    x=x0+q(:,1:k)*y;
    r=norm(b-A*x);
    [k r]
end
x
