## 4.6-4
##
source('utils.R')

A = matrix(c(7, 8, 9, 8, 9, 10, 9, 10, 8), nrow=3, ncol=3, byrow=T)
print(A)

B = matrix(c(24, 27, 27, 24.1, 26.9, 26.9), nrow=3, ncol=2, byrow=FALSE) ## the two right hand sides 
print(B)

## Solve A X = B
##
res = GE_w_scaled_partial_pivoting(A, B)
X1 = res$X
print(X1)

## The updated solution due to iterative improvement:
##
R = B - A %*% X1

## Solve A E = R for E
##
res = GE_w_scaled_partial_pivoting(A, R)
E = res$X

X2 = X1 + E
print(X2)

## Estimate the condition number:
##
x1_minus_x2 = X2[, 1] - X2[, 2]
x1_minus_x2 = matrix(data=x1_minus_x2, nrow=3, ncol=1, byrow=T)

b1_minus_b2 = B[, 1] - B[, 2]
b1_minus_b2 = matrix(data=b1_minus_b2, nrow=3, ncol=1, byrow=T)

my_kappa = norm(x1_minus_x2, '2')/norm(b1_minus_b2, '2')
print(sprintf('estimated kappa= %f', my_kappa))

## What is the "true" condition number for A:
##
print(kappa(A, norm='2'))