%this routine performs qr factorization using householder matrices
%it is correct, but highly inefficient. In assignment 4, cs513 students
%will examine its inefficiency and will improve it.

%the file is modified version of sloppy_qr.m: many pause, disp
% were added, and out was not suppressed any more.

disp(' n=size(A)              % n is defined as the order of A; ')
  n=size(A)              % n is defined as the order of A;
  pause
disp(' R=A;                    ')
  R=A;                    

disp(' Q=eye(n(1)); ')
  Q=eye(n(1));
disp(' for i=1:min(n(1)-1,n(2)), ')
  for i=1:min(n(1)-1,n(2)),

%i is a running index for the elimination step. entering
%the i'th step, we have in hand a factroization
%Q*R, with Q the product of the previous Householder
%matrices we used, and R the partially processed matrix (i.e., A_i in the
%notations used in class). During the ith step,
%we generate a new householder matrix H and rewrite Q*R as
%(Q*H)*(H*R). We then replace R H*R and Q by Q*H.
%Initially R is the given matrix, and Q is the identity matrix
disp(' x=R(:,i)                        %  x is the current i th column of R ')
     x=R(:,i)                          %  x is the current i'th column of R
disp(' a=norm(x(i:n),2);     ')
     a=norm(x(i:n),2);     
disp(' y=[x(1:i-1)'' a zeros(1,n-i)]''  % y is the desired i th column of R ')
     y=[x(1:i-1)' a zeros(1,n-i)]'  % `y' is the desired i'th column of R
     pause
disp(' w=(x-y)/norm(x-y,2)            % w is the vector that determines H ')
     w=(x-y)/norm(x-y,2)            % w is the vector that determines H
     pause
     !C
disp(' H=eye(n(1))-2*w*w'' ')
     H=eye(n(1))-2*w*w'
     pause
     !C
disp(' Q=Q*H, R=H*R ')
     Q=Q*H, R=H*R
     pause
  end
!C
disp(' A, Q, R ')
A, Q, R
pause
Q=Q(:,[1:n(2)]), R=R([1:n(2)],:)
pause
disp(' norm(A-Q*R), norm(eye(n(1))-Q''*Q) ')
  norm(A-Q*R), norm(eye(n(2))-Q'*Q)

