## 4.4-2
##
source('utils.R')

A = matrix(data=c(0.21, 0.32, 0.12, 0.31, 0.1, 0.15, 0.24, 0.22, 0.2, 0.24, 0.46, 0.36, 0.61, 0.4, 0.32, 0.2), nrow=4, ncol=4, byrow=TRUE)
print('A=')
print(A)

## Compute the inverse of A by computing Gaussian Elimination on A x = e_j for 1 <= j <= n:
##
n = dim(A)[1]
AInv = matrix(NA, nrow=n, ncol=n, byrow=TRUE)
for( jj in 1:n ){
    b = matrix(0, nrow=n, ncol=1)
    b[jj] = 1
    results = GE_w_scaled_partial_pivoting(A, b)
    AInv[, jj] = results$X
}
print('AInv=')
print(AInv)

print('A times AInverse:')
print(A %*% AInv)
print('AInverse times A:')
print(AInv %*% A)

## Compute the inverse of A by doing a'factorization'and then several'back-substitutions'
##
results = GE_w_scaled_partial_pivoting(A, diag(n))
AlInv2 = results$X
print(AInv - AlInv2)