% Skripta testira razne algoritme za QR-faktorizaciju.
% Generira 100 matrica visoke kondicije (cnd) i za svaku pokrene CGS, MGS, VVt i ugradjeni MATLABov qr.
% Za dobivene faktorizacije A = QR nacrta grafove na kojima je prikazano ||A-Q*R||, odnosno ||Q'*Q - I||.

numEx = 100;

cgsRez = zeros( 100, 1 );
mgsRez = zeros( 100, 1 );
VVtRez = zeros( 100, 1 );
qrRez  = zeros( 100, 1 );

cgsRez2 = zeros( 100, 1 );
mgsRez2 = zeros( 100, 1 );
VVtRez2 = zeros( 100, 1 );
qrRez2  = zeros( 100, 1 );


for ex = 1 : numEx
	A = randn(100);
	[u, s, v]=svd(A);
	% Promijenimo singularne vrijednosti tako da variraju od 1 to cnd.
	cnd = 100000;
	sd = [1, cnd, exp(rand(1,100-2)*log(cnd))];
	s = diag(sd);
	A = u * s * v';
	
	[Q, R] = cgs( A );
	cgsRez( ex ) = norm( Q'*Q - eye( 100 ) );
	cgsRez2( ex ) = norm( Q*R - A );
	
	[Q, R] = mgs( A );
	mgsRez( ex ) = norm( Q'*Q - eye( 100 ) );
	mgsRez2( ex ) = norm( Q*R - A );
	
	[Q, R] = VVt( A );
	VVtRez( ex ) = norm( Q'*Q - eye( 100 ) );
	VVtRez2( ex ) = norm( Q*R - A );
	
	[Q, R] = qr( A, 0 );
	qrRez( ex ) = norm( Q'*Q - eye( 100 ) );
	qrRez2( ex ) = norm( Q*R - A );
end

semilogy( cgsRez, 'r' ); hold on;
semilogy( mgsRez, 'b' ); 
semilogy( VVtRez, 'g' ); 
semilogy( qrRez, 'k' );

legend( 'CGS', 'MGS', 'q=a-Q*Qt*a', 'MATLAB QR' );
title( sprintf( 'norm(Q*Qt - I) za matrice A kondicije %d', cnd ) );
hold off;

figure;
semilogy( cgsRez2, 'r' ); hold on;
semilogy( mgsRez2, 'b' ); 
semilogy( VVtRez2, 'g' ); 
semilogy( qrRez2, 'k' );

legend( 'CGS', 'MGS', 'q=a-Q*Qt*a', 'MATLAB QR' );
title( sprintf( 'norm(Q*R - A) za matrice A kondicije %d', cnd) );
hold off;
