function ex03()
%% a.
r = @(x) runde(x, 4);
A = [50 -20 91; -29 -59 10; 110 90 81];
P1 = [0 0 1; 0 1 0; 1 0 0];
A1 = P1*A

q1 = r(A1(2,1)/A1(1,1))
z2 = r(A1(2,:) - r(q1*A1(1,:)));
q2 = r(A1(3,1)/A1(1,1))
z3 = r(A1(3,:) - r(q2*A1(1,:)));
A2 = [A1(1,:); z2; z3]

P2 = [1 0 0; 0 0 1; 0 1 0];
A3 = P2*A2
q3 = r(A3(3,2)/A3(2,2))
z3 = r(A3(3,:) - r(q3*A3(2,:)));
U = [A3(1,:); A3(2,:); z3]

L = [1 0 0; q2 1 0; q1 q3 1]
P = P2*P1
P*A
r(L*U)

%% b.
x = [4 1 -1]';
b = [89 -185 449]';

c = P*b

y1 = c(1);
y2 = r(c(2) - r(L(2,1)*y1));
y3 = r(c(3) - r(L(3,1)*y1) - r(L(3,2)*y2));
y = [y1 y2 y3]'

x3 = r(y(3)/U(3,3));
x2 = r(r(y(2) - r(U(2,3)*x3))/U(2,2));
x1 = r(r(y(1) - r(U(1,2)*x2) - r(U(1,3)*x3))/U(1,1));
xS = [x1 x2 x3]'

rB = norm(A*xS - b)
eB = norm(x - xS)/norm(x)

%% c.
[L, U, P] = lu(A)
y = L \ (P*b)
xSc = U \ y

rC = norm(A*xSc - b)
eC = norm(x - xSc)/norm(x)
condA = cond(A)

rhsC = cond(A)*rC/norm(b)
rhsB = cond(A)*rB/norm(b)
