function Gram_Schmidt (b,flag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Gram_Schmidt.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%%% demonstration of the classical and modified 
%%% Gram-Schmidt orthogonalization algorithms 
%%% refs: [Trefethen, Bau, "Numerical Linear Algebra"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2004/10
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% INSTRUCTIONS: this is a function that can be called 
%% input: matrix b 
%% output: the algorithms produce: 
%% an upper triagular array r of factors
%% and an array q of orthogonal columns 
%% if flag ==0 only the orthogonality is checked
%% if flag !=0 also the columns of QT*Q are shown
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if nargin == 0 return; end;
if nargin == 1 flag = 0; end; 

m = size(b,1); 
n = size(b,2); 

q = zeros(m,n); r = zeros(n);
%%% classical Gram-Schmidt
for j = 1:n
  v(:,j) = b(:,j);
  for i=1:j-1
    r(i,j) = q(:,i)' * b(:,j);
    v(:,j) = v(:,j) - r(i,j) * q(:,i);
  end;
  r(j,j) = norm(v(:,j));
  q(:,j) = v(:,j) / r(j,j);
end;
class_r = r; class_q = q;

q = zeros(m,n); r = zeros(n);
%%% modified  Gram-Schmidt
for i = 1:n
  v(:,i) = b(:,i);
end;
for i = 1:n
  r(i,i) = norm(v(:,i));
  q(:,i) = v(:,i) / r(i,i);
  for j = i+1:n
    r(i,j) = q(:,i)'*v(:,j);
    v(:,j) = v(:,j) - r(i,j)*q(:,i);
  end;
end;
mod_r = r; mod_q = q;

format long;
fprintf('Results of classical Gram=Schmidt: loss of ortogonality ||QT*Q||=%g\n',...
        norm (class_q'*class_q-eye(n)));
if (flag>0) class_q'*class_q,pause; else pcolor(class_q);title('Classical Q');pause,end;

fprintf('Results of modified Gram=Schmidt:  loss of ortogonality ||QT*Q||=%g\n',...
        norm (mod_q'*mod_q-eye(n)));
if (flag>0) mod_q'*mod_q,pause;     else pcolor(mod_q);title('Modified Q');pause,end;
format;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
