% Demonstration of QR factorization with Gram-Schmidt 
% orthogonalization.  H. Motteler, Feb 1999
%
% This example is intended for transparency, rather than
% efficiency.  The semicolon at the end of a statement can
% be deleted to watch the value change with each iteration.


% The array A and vector b from Example 3.5 of the text
%
A0 = [1   -1    1;
     1   -0.5  0.25;
     1    0    0;
     1    0.5  0.25;
     1    1    1];

b0 = [1 0.5, 0, .5, 2]';


% a larger array A and vector b, for a more serious teset
%
% A0 = rand(20, 18);
% b0 = rand(20, 1);


% In the following modified Gram-Schmidt orthogonalization,
% the matrix A is changed in place to become Q, or actually
% Q1 of a partitioned QR factorization.  

A = A0
b = b0;
[m, n] = size(A);
R = zeros(n,n);

for k = 1 : n

  R(k,k) = norm(A(:,k))

  A(:,k) = A(:,k) / R(k,k)

  for j = k+1 : n

    R(k,j) = A(:,k)' * A(:,j);

    A(:,j) = A(:,j) - R(k,j) * A(:,k);

  end
end

% check residual
%
Q1 = A;
norm(A0 - Q1*R)

% to solve A*X = b for x, we can now solve the upper
% triangular system R*x = Q'*b by back substitution