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

A = matrix(data=c(2, 2, 1, 1, 1, 1, 3, 2, 1), nrow=3, ncol=3, 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)

print('Needed permuation of the rows for a LU factorization ps=')
print(results$ps)
print('PA=')
PA = A[results$ps, ]
print(PA)

results = GE_w_scaled_partial_pivoting(PA, diag(n)) ## this call should need no permutations so W = [L, U]
print(results)
print('Product LU=')
print(results$L %*% results$U)