function testgram (alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% testgram.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%%% demonstration of how the classical and modified 
%%% Gram-Schmidt orthogonalization algorithms 
%%% using examples by Lauchle and Bjorck
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: call as testgram
%% or testgram(alpha) where alpha should be a small number
%% uses Gram_Schmidt.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska, 2004/11-2005/03
%% D. Fettig  2005/03
%% supported by DeLoach Work Scholarship Program 
%% of the Oregon State University Honors College
%% http://oregonstate.edu/dept/honors/students/opportunities/DeLoach.htm
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if nargin == 0 alpha = sqrt(eps); end;

fprintf('Demonstration of the loss of orthogonality by classical Gram-Schmidt        ');
fprintf('and its recovery by the modified algorithm using the example by Lauchli 1959');
%
b = zeros(4,3);
b (:,1) = [1, alpha, 0, 0]';
b (:,2) = [1, 0, alpha, 0]';
b (:,3) = [1, 0, 0, alpha]';
%
b
%
Gram_Schmidt (b,1);

%%%% this function can be tested with examples from paper by [Rice, 1966]
clear b;
fprintf('Examples from paper by Rice, 1966\n');
%
p = inline('x.^3+x.^2+x+1');
n = 20; m = 15;
for i=1:n for j=1:m b(i,j) = p(j) * sin(i*j*pi/(2*n));end;end;
pcolor (b),title('Original matrix');pause,
Gram_Schmidt (b);
