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

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

%**** COMPUTER PROBLEM 2.1 *****
%(a) Show that the matrix A is singular

A=[0.1 0.2 0.3; 0.4 0.5 0.6; 0.7 0.8 0.9];
R=rank(A)
AI=inv(A)

% Describe the set of solutions for:

b=[0.1 0.3 0.5]';
x=A\b

%Answer (a): There could be no solution for b or infinite many solutions.

%(b) At what point would the  Gaussian elimination with partial pivoting 
% in exact arithmetic fail?

echo off;
n=3;
M=eye(n);	% The place for the matrix M (start with the identity matrix),
A1=A;		% and original matrix.
j=1;		% first pass	
M1=eye(n);
P1=[0 0 1;0 1 0;1 0 0]
A1=P1*A1	% interchange 3. and 1. row 
   for i=j+1:n
      M1(i,j)=-A1(i,j)/A1(j,j);
   end
M=M1*M
A1=M1*A1
j=2;		% second pass	
M2=eye(n);
P2=[1 0 0;0 0 1;0 1 0]
A2=P2*A1	% interchange 3. and 2. row 
   for i=j+1:n
      M2(i,j)=-A2(i,j)/A2(j,j);
   end
M=M2*M
A2=M2*A2
echo on;


%Answer (b): The Gaussian elimination would not fail, because the zero pivot  
%            takes place still in the final third pass.

%(c) Calculate cond(A). Estimate the number of accurate digits 
%in the solution.

cond(A)
log10(cond(A))

% Answer (c):
% We can expect the loss of 16 digits in the result. The result is arbitrary.


%**** COMPUTER PROBLEM 2.15 *****
% Implement Gaussian elimination using each of the six different orderings
% of triple nested loop, and compare performances.

%%% generic implementation
n=4;
A=rand(n,n)
b=rand(n,1)
echo off;
M=eye(n);	% The place for the matrix M (start with the identity matrix),
A1=A;		% and original matrix.
for j=1:n-1
   M1=eye(n);
   for i=j+1:n
      M1(i,j)=-A1(i,j)/A1(j,j);
   end
   M=M1*M;
   A1=M1*A1;
end
echo on;
M		% Gaussian elimination matrix M.
U=A1		% Modified matrix A is upper triangular.
L=inv(M)	% Inverse of M is lower triangular.
L*U		% Should be the same as original A.
x=(A\b)	% The  MatLab solution of the original random system.
bm=(M*b)	% Generate the modified right-hand vector and solve modified system.
x=(A1\bm)	% The solution of modified system should be the same as MatLab solution.

%%%Usually, the Gaussian elimination can be implemented as a triple loop.
% The matrix M and L are not generated explicitly. The transformed matrix (upper triangular)
% is saved in the upper triangle of the original matrix A, the transformed right-hand vector
% is saved on the place of the original right-hand vector b.

echo off;
bg=b;		% Make a copy of b and A
Ag=A;
for k=1:n-1
   for i=k+1:n
      M(i,k)=Ag(i,k)/Ag(k,k);
      for j=k+1:n
         Ag(i,j)=Ag(i,j)-M(i,k)*Ag(k,j);
      end
      bg(i)=bg(i)-M(i,k)*bg(k);
   end
end
echo on;
Ag		% Compare upper triangle of Ag with U.
bg		% Compare bg with bm.

%%% replace indices k and i. The result should be the same.
echo off;
bg=b;		% Make a copy of b and A
Ag=A;
for i=2:n		% replace and change range
   for k=1:i-1		% replace and change range
      % M(i,k)=Ag(i,k)/Ag(k,k); % bring in inner loop
      for j=k+1:n  	
         Ag(i,j)=Ag(i,j)-(Ag(i,k)/Ag(k,k))*Ag(k,j);
      end
      % bg(i)=bg(i)-M(i,k)*bg(k);
   end
end
echo on;
Ag		% Compare upper triangle of Ag with U.


%**** COMPUTER PROBLEM 2.16 *****
% Implement forward and backward substitution for different order of
% two indices in nested loops(2 variants for each algorithm).
% Use the same system as in Computer problem 2.15.

%%%%%Forward substitution for the lower triangular matrix
L
echo off;
xf=zeros(n,1);		% Initialize the solution vector
xf(1)=b(1)/L(1,1);	%First value
for i=2:n 
   for j=1:i-1
      xf(i)= xf(i) + L(i,j)*xf(j);
   end
   xf(i) = (b(i)-xf(i))/L(i,i);
end
echo on
xf		% The forward substitution (standard index order).
y=L\b		% The solution of the upper triangular system should be the same 
		% as by forward substitution.


%%%%%Backward substitution for the upper triangular matrix
U
echo off;
xb=zeros(n,1);		% Initialize the solution vector
xb(n)=b(n)/U(n,n);	%First value
for i=(n-1):-1:1 
   for j=(i+1):n
      xb(i)= xb(i) + U(i,j)*xb(j);
   end
   xb(i) = (b(i)-xb(i))/U(i,i);
end
echo on
xb		% The backward substitution (standard index order).
xo=U\b		% The solution of the upper triangular system should be the same
		% as by backward substitution.
x=U\y		% The solution of the complete system should be the same
		% as in the Computer problem 2.15.


%**** UPGRADE OF COMPUTER PROBLEM 2.16 *****
% Upgrade computer problem 2.16 with the following options:
% From standard Gaussian elimination derive a routine for solving tridiagonal 
% system. Then include partial pivoting. Adapt the program also for positive 
% symmetric positive definite systems and write Cholesky factorization.

n=5;			% System order
bc=rand(n,1)		% right-hand-side vector
x=bc;
d=rand(n,1);		% main diagonal
du=rand(n-1,1);		% upper diagonal
dl=rand(n-1,1);		% lower diagonal
Ac=diag(d)+diag(du,1)+diag(dl,-1)	% Tridiagonal matrix
X=Ac;

%%%% Gaussian elimination for tridiagonal system
echo off;
for k=1:n-1		% for all columns
   for i=k+1:k+1	% for all remaining rows in the current column (only one nonzero row)
      M(i,k)=Ac(i,k)/Ac(k,k);	% calculate current column of the elementary elimination matrix
      for j=k+1:k+1	% annihilate remaining entries in current column (only a single entry)
         Ac(i,j)=Ac(i,j)-M(i,k)*Ac(k,j);
      end
      bc(i)=bc(i)-M(i,k)*bc(k);
   end
end
echo on;
% We can discard two inner loops for tridiagonal systems because there remains only a 
% single pass i=j=k+1. The calculation complexity remains n.
U=Ac
b=bc

%%%% Backward substitution for the upper triangular part of Ac gives the solution.
echo off;
xb=zeros(n,1);		% Initialize the solution vector
xb(n)=b(n)/U(n,n);	%First value
for i=(n-1):-1:1		% for all columns 
   for j=(i+1):(i+1)	% only for one additional element
      xb(i)= xb(i) + U(i,j)*xb(j);
   end
   xb(i) = (b(i)-xb(i))/U(i,i);
end
echo on
xb		% The solution by backward substitution should be the same as 
xo=X\x		% the solution of the original problem. The calculation complexity remains n.


% If Pivoting is necessary, then the largest element in the column (b or a) has to be chosen
% for pivot.

% If both side diagonal are the same and the matrix is positive definite the Cholesky 
% factorization can be implemented.
d=[1:n];		% make the matrix positive definite by the main diagonal
du=rand(n-1,1);		% upper diagonal
Ac=diag(d)+diag(du,-1)	% Take only the lower part.

% Implement Cholesky factorization
echo off;
for j=1:n		% for each column j
   for k=j-1:j-1 if (k~=0)	% only for one prior column k!
   		% In general for all prior columns: for k=1:j-1	      
      for i=j:j+1	if (i<=n) % for two elements from the main and low diagonal of column j
         			% in general for all elements of column j below diagonal: for i=j:n
          Ac(i,j)=Ac(i,j)-Ac(i,k)*Ac(j,k); end
      end; end
   end
   Ac(j,j)=sqrt(Ac(j,j));
   for k=j+1:j+1 if (k<=n)	% Scale one additional element from the lower diagonal.
      				% In general all remaining elements below diagonal: for k=j+1:n
       Ac(k,j)=Ac(k,j)/Ac(j,j); end
   end
end
echo on;
Ac		% Ac is equal to L, the overall complexity n*2+n*2=4*n
Ac*Ac'		% Yields the original matrix.

% To get the solution x would proceed by forward and backward substitution of LL'x=b
% for two systems: Ly=b and L'x=y.