% Numerik D-MAVT, RW (Nipp), SS 06
%
% This is the M-file solution to the second exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
% Then to run the file type "s2" in the command window
%
function s2
% Serie 2, A1a)
disp(['Aufgabe 1a']);
format compact;
A=[.005 1; 1 1]; y=[.5;1];
x=A\y,
% output:
% x =
%    0.50251256281407
%    0.49748743718593
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Serie 2, A1b)
disp(['Aufgabe 1b']);
format compact;
A=[1 1 2;-1 2 1;0 -1 1]; b=[1;2;3];
[L,U,P]=lu(A), % P*A=L*U, if P==identity: no permutations made!
b=P*b;	% permute rhs (actually not necessary here)
y=L\b,
x=U\y,
% output:
% L =
%     1.0000         0         0
%    -1.0000    1.0000         0
%          0   -0.3333    1.0000
% U =
%      1     1     2
%      0     3     3
%      0     0     2
% P =
%      1     0     0
%      0     1     0
%      0     0     1
% y =
%      1
%      3
%      4
% x =
%     -2
%     -1
%      2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function s3
% Serie 2, A2
disp(['Aufgabe 2']);
format long; 
% to check convergence look at spectral radius of iteration matrix T
A=[2 1 0; 1 2 1; 0 1 2]; 			% the matrices to test
B=[1 -1 0; 5 2 0; 0 0 1];
C=[3 2 1; 2 3 2; 1 2 3];
DA=diag(diag(A)); LA=tril(A)-DA; RA=triu(A)-DA;      % A=DA+LA+RA
DB=diag(diag(B)); LB=tril(B)-DB; RB=triu(B)-DB;      % B=DB+LB+RB
DC=diag(diag(C)); LC=tril(C)-DC; RC=triu(C)-DC;      % C=DC+LC+RC
TAj=-inv(DA)*(LA+RA); 	% x^{k+1}=T*x^k+c jacobi 
TBj=-inv(DB)*(LB+RB); 	
TCj=-inv(DC)*(LC+RC); 
% if rho < 1, we have convergence (assume (I-T) is regular)
rhoAj=max(abs(eig(TAj))),
rhoBj=max(abs(eig(TBj))),
rhoCj=max(abs(eig(TCj))),
% output:
% rhoAj =
%    0.70710678118655
% rhoBj =
%    1.58113883008419
% rhoCj =
%    1.12409377442300
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Serie 2, A3
% test with A1a
disp(['Aufgabe 3, testing with A1a']);
A=[0.005 1; 1 1]; [L,R,P]=lrp(A),
% output:
% L =
%     1.0000         0
%     0.0050    1.0000
% R =
%     1.0000    1.0000
%          0    0.9950
% P =
%      0     1
%      1     0

% test with A1b
disp(['Aufgabe 3, testing with A1b']);
A=[1 1 2;-1 2 1;0 -1 1]; [L,R,P]=lrp(A),
% output:
% L =
%     1.0000         0         0
%    -1.0000    1.0000         0
%          0   -0.3333    1.0000
% R =
%      1     1     2
%      0     3     3
%      0     0     2
% P =
%      1     0     0
%      0     1     0
%      0     0     1
return;



% Serie 2, A3 
function [L,R,P] = lrp(A)

% desc:	[L,R,P] = lrp(A) stores an upper triangular matrix in R and a
%	lower triangular matrix in L, so that P*A = L*R.  A must be square.
%	relative column maximal pivoting, error if any pivot is 0
% auth: fab@math.ethz.ch

n=size(A); n=n(1);	% A is n \times n
L=eye(n);  R=zeros(n); P=eye(n);
for k=1:n-1 		% loop over rows (without last)
  % pivoting:
  m=max(abs(A(k:n,k:n))')';	% maximal element of each row 
  q=abs(A(k:n,k))./m;		% quotient
  [qp,p]=max(q);		
  if (qp==0) error('only zero pivots'); end;		
  p=p+k-1;			% take a_{pk} as pivot
  if (p~=k)			% we have to swap rows in A and P
    disp(['k=',int2str(k),': swap p<->k: ',int2str(p),'<->',int2str(k)]);
    v=A(k,:); w=P(k,:);		% swap rows p <--> k
    A(k,:)=A(p,:); P(k,:)=P(p,:); 
    A(p,:)=v;  P(p,:)=w;  
    % swap old L_ik's
    u=L(k,1:k-1); L(k,1:k-1)=L(p,1:k-1); L(p,1:k-1)=u;
  end
  % end of pivoting
  for i=k+1:n 		% loop over columns
    L(i,k) = A(i,k)/A(k,k);
    A(i,k:n) = A(i,k:n) - A(k,k:n)*L(i,k);	% get new row
  end
end
R=A;
return;
