% Test script that shows how Modified Gram-Schmidt (MGS) can be viewed as a
% sequence of triangular transformations applied from the right of the
% matrix.  
%
% You are encouraged to play around with the matrices R1, R2, R3, R4.
% What is inv(R1)?
% What is inv(R2)?
% What is inv(R2)*inv(R1)?

function main_mgs_triangular

rng(0)
A = randn(4);

r11     = norm(A(:,1));
q1      = (1/r11)*A(:,1);
r12     = q1'*A(:,2);
r13     = q1'*A(:,3);
r14     = q1'*A(:,4);
R1      = eye(4);
R1(1,:) = (1/r11)*[1,-r12,-r13,-r14];
A1      = A*R1;

keyboard

r22     = norm(A1(:,2));
q2      = (1/r22)*A1(:,2);
r23     = q2'*A1(:,3);
r24     = q2'*A1(:,4);
R2      = eye(4);
R2(2,:) = (1/r22)*[0,1,-r23,-r24];
A2      = A1*R2;

keyboard

r33     = norm(A2(:,3));
q3      = (1/r33)*A2(:,3);
r34     = q3'*A2(:,4);
R3      = eye(4);
R3(3,:) = (1/r33)*[0,0,1,-r34];
A3      = A2*R3;

keyboard

r44     = norm(A3(:,4));
R4      = eye(4);
R4(4,4) = (1/r44);
A4      = A3*R4;

Q       = A4;
R       = inv(R4)*inv(R3)*inv(R2)*inv(R1);

fprintf(1,'|A - Q*R|        = %10.3e\n',max(max(abs(A - Q*R))))
fprintf(1,'|I - trans(Q)*Q| = %10.3e\n',max(max(abs(eye(4) - Q'*Q))))

keyboard


return