## 4.2-8
##
## Work this problem "by hand" i.e. perform Algorithm 4.2 "by hand"
##

n_digits = 3

## Initialize the W matrix:
##
W = matrix(data=c(0.21, 0.32, 0.12, 0.31, 0.96, 0.1, 0.15, 0.24, 0.22, 0.71, 0.2, 0.24, 0.46, 0.36, 1.26, 0.61, 0.4, 0.32, 0.2, 1.53), nrow=4, ncol=5, byrow=TRUE)
W_orig = W
print('Initial W matrix:')
print(W_orig)

n = dim(W)[1]

kk = 1

## The second row:
##
ii = 2
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print(W)

## The third row:
##
ii = 3
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print(W)

## The fourth row:
##
ii = 4
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print(W)

kk = 2

## The third row:
##
ii = 3
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print(W)

## The fourth row:
##
ii = 4
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print(W)

kk = 3

## The fourth row:
##
ii = 4
m = W[ii, kk] = signif(W[ii, kk] / W[kk, kk], n_digits)
jj = (kk+1):(n+1)
W[ii, jj] = signif(W[ii, jj] - m * W[kk, jj], n_digits)
print('final W before backsubstitution:')
print(W)

x = matrix(data=0, nrow=n, ncol=1, byrow=TRUE)

## Backsubstitution:
##
kk = 4
x[kk] = signif(W[kk, n+1]/W[kk, kk], n_digits)
print(x)

kk = 3
x[kk] = signif((W[kk, n+1] - sum(W[kk, (kk+1):n] * x[(kk+1):n]) )/W[kk, kk], n_digits)
print(x)

kk = 2
x[kk] = signif(( W[kk, n+1] - sum(W[kk, (kk+1):n] * x[(kk+1):n]) )/W[kk, kk], n_digits)
print(x)

kk = 1
x[kk] = signif((W[kk, n+1] - sum(W[kk,(kk+1):n] * x[(kk+1):n]))/W[kk, kk], n_digits)
print(sprintf('n_digits in signif= %d', n_digits))
print(x)