echo on
%%%%% ch03_pro.m %%%%%
% Run script and study the effects of each line.

format compact
%************************************************
%**** Some Computer Problems from Chapter 03 ****
%************************************************

%**** COMPUTER PROBLEM 3.x *****
% Factorize a selected square matrix by:
%(a) LU by Gaussian elimination  without pivoting
%(b) LL' by Cholesky factorization if the matrix is symmetric and positive definite
%(c) QR by Householder transformation
%(d) QR by Givens Rotation and
%(e) QR by Gram-Schmidt orthogonalization.
% Compare the calculation complexity for all above methods.

%%%% Input data
%t=[-1.0 -0.5 0.0 0.5 1.0]; A=[[ones(1,5)]' t' t'.^2] % Five data points and Vandermonde matrix for quadratic polynomial (5x3)
%n=4; A=pascal(n) 
n=4; A=rand(n,n)
[m n]=size(A);	


%%%% Implement Gaussian LU factorization without pivoting (only for square matrices ???)
echo off;
Ag=A(1:n,:); % Take square matrix for input
Ga=eye(n);	% Gaussian transformation matrix usually not needed explicitly
flops(0);
for k=1:n-1	% for each column k
   for i=k+1:n	% scale all remaining elements of column k below diagonal by pivot Ag(k,k)  
      Ag(i,k)=Ag(i,k)/Ag(k,k);	% and leave element of L on position below diagonal 
      				% this elements would be annihilated in U, but belong to L
      for j=k+1:n	% for each scaled element, update all row i elements of former A,
         Ag(i,j)=Ag(i,j)-Ag(i,k)*Ag(k,j);  % ie. future U. 
      end
   end
end
echo on;
Ag	% Upper triangle of Ag is U and lower subdiagonal entries plus I is L. 
	% The system Ax=b transforms, because L*U=A, in the form L*U*x=b.
flops	% Cost of LU factorization is O(n^3/3), no additional memory needed.


%%%% Implement Cholesky LL' factorization
echo off;
Ac=A'*A;		% Make a symmetric positive definite matrix A'*A
flops(0);
for j=1:n		% for each column j
   for k=1:j-1 		% for all prior columns k	      
      for i=j:n		% update all elements of column j below diagonal
          Ac(i,j)=Ac(i,j)-Ac(i,k)*Ac(j,k); 
      end
   end
   Ac(j,j)=sqrt(Ac(j,j));
   for k=j+1:n	 	% scale all remaining elements below diagonal
       Ac(k,j)=Ac(k,j)/Ac(j,j); % lower triangle equals L
       Ac(j,k)=Ac(k,j);	% Upper triangle equals L' 
   end
end
echo on;
Ac	% Lower triangle of Ac is equal to L of A, upper triangle is L'
	% so that L*L'=A. The system: A*x=b yields two triangular systems L*y=b, L'*x=y.
flops	% Cost of LU factorization is O(1/2*n^3/3), no additional memory needed.


%%%% Implement Householder QR factorization
echo off;
Ah=A;		% Make a copy of A
Hh=eye(m);
flops(0);
for j=1:n	% for all columns j form a vector v 
   v=Ah(:,j);	% make a copy of column j
   for k=1:j-1 v(k)=0;	% all elements in v above diagonal should be 0 
   end
   no=norm(v);	% calculate norm of vector v
	if (v(j)<0) v(j)=v(j)-no; % take the same sign as A(j,j) to prevent cancellation
   else v(j)= v(j)+no; 
   end
   H=eye(m)-(2*(v*v')/(v'*v)); % implement elementary Householder transformation if needed
   Hh=H*Hh;		% Build the optional transformation matrix Q' by premultiplying.
   for i=j:n	% for all unfactorized columns i
      Ah(:,i)=Ah(:,i)-(2*(v'*Ah(:,i))/(v'*v))*v; % implement Householder transformation
   end
end
echo on;
flops	% Cost of Householder QR factorization is O(5*n^3/3)
Hh	% Householder transformation matrix Hh=Q'(in general not needed for orthogonalization)
Ah	% Upper triangle of Ah is R. 
Hh*A	% Q'*A=R and system Ax=b becomes Q*R*x=b (solution: Qy=b and Rx=y)
Hh'*Hh	% Hh is orthogonal transform., hence Q*Q'=I so A=Q*R. System Ax=b becomes R*x=Q'*b.
norm(A)	% Householder transformation preserves norm and hence the solution of the system.
norm(Ah)  


%%%% Implement Givens rotation for QR factorization 
echo off;
Ai=A;		% Make a copy of A
G=eye(m,m);	% Place for givens transformation
flops(0);
for j=1:n	% for all columns j annihilate 
   for k=m:-1:j+1	% subdiagonal entries k starting from bottom
      sq=sqrt(Ai(k,j)^2+Ai(k-1,j)^2);	
      c=Ai(k-1,j)/sq;	% calculate rotation c
      s=Ai(k,j)/sq;	% and s
      Gi=eye(m);
      Gi(k-1,k-1)=c; Gi(k,k)=c; %embedded rotations in elementary transformation
      Gi(k,k-1)= -s; Gi(k-1,k)=s;
      G=Gi*G;	% Build transformation matrix
      Ai=Gi*Ai;	% Build upper triangular modified A (modified right-hand side Gi*b, if needed)
   end
end
echo on;
Ai	% Upper triangle of Ai=G*A is equal to R. Q is R' and Q*R=A.
flops	% Cost of Givens rotation for QR factorization is O(5*n^3/3)
G'*G	% Givens transformation is orthogonal, 
norm(A)	% so it preserves norm and the solution of the system.
norm(Ai)  


%%%% Implement modified Grahm-Schmidt orthogonalization
echo off;
As=A;		% Make a copy of A, for future place for Grahm-Schmidt transformation
Rs=zeros(n,n);	% prepare extra space for upper triangular R
flops(0);
for k=1:n	% for each column k
   Rs(k,k)=norm(As(:,k));     % calculate norm
   As(:,k)=As(:,k)/Rs(k,k);   % normalized AS(:,k) overwrites original AS(:,k)
   for j=k+1:n
      Rs(k,j)=As(:,k)'*As(:,j);   % orthogonalize against all subsequent columns
      As(:,j)=As(:,j)-Rs(k,j)*As(:,k);   % and update them
   end
end
	   
echo on;
flops	% Cost of Grahm-Schmidt rotation for QR factorization is O(5*n^3/3)
As  % Gram-Schmidt transformation Q
Rs	% Rs is equal to R. Q*R=A so Ax=b can be solved as R*x=Q'*b
As'*As	% Grahm-Schmidt transformation is orthogonal, 
norm(A)	% so it preserves norm and the solution of the system.
norm(As*A)	% and normalized