## 4.3-3:
##

n_digits = 4

x_true = c(1, 0.25)

if( T ){
    ##
    ## Part (a): Using the first equation as a pivot (i.e. perform Algorithm 4.2 "by hand" using the original row ordering):
    ##
    W = matrix(data=c(0.141e-2, 0.4004e-1, 0.1142e-1, 0.2, 0.4912e1, 0.1428e1), nrow=2, ncol=3, byrow=TRUE) ## Initialize the W matrix
    print('Initial W matrix (using the first equation as a pivot):')
    print(W)

    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('Final W before backsubstitution:')
    ##print(W)

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

    ## Backsubstitution:
    ##
    kk = 2
    x_a[kk] = signif(W[kk, n+1]/W[kk, kk], n_digits)

    kk = 1
    x_a[kk] = signif((W[kk, n+1] - sum(W[kk,(kk+1):n] * x_a[(kk+1):n]))/W[kk, kk], n_digits)
    print(sprintf('Using the first equation as a pivot (n_digits in signif= %d); x_a=', n_digits))
    print(x_a)
    print(sprintf('error in x_a= %10.6f', sum((x_a-x_true)^2)))
}


if( T ){
    ##
    ## Part (b): Using the second equation as a pivot (i.e. perform Algorithm 4.2 "by hand" after exchanging the original row ordering):
    ##
    W = matrix(data=c(0.2, 0.4912e1, 0.1428e1, 0.141e-2, 0.4004e-1, 0.1142e-1), nrow=2, ncol=3, byrow=TRUE) ## Initialize the W matrix
    print('Initial W matrix (using the second equation as a pivot):')
    print(W)

    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('Final W before backsubstitution:')
    ##print(W)

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

    ## Backsubstitution:
    ##
    kk = 2
    x_b[kk] = signif(W[kk, n+1]/W[kk, kk], n_digits)

    kk = 1
    x_b[kk] = signif((W[kk, n+1] - sum(W[kk,(kk+1):n] * x_b[(kk+1):n]))/W[kk, kk], n_digits)
    print(sprintf('Using the second equation as a pivot (n_digits in signif= %d); x_b=', n_digits))
    print(x_b)
    print(sprintf('error in x_b= %10.6f', sum((x_b-x_true)^2)))
}

if( T ){
    ##
    ## Part (c): Performing total pivoting:
    ##
    W = matrix(data=c(0.141e-2, 0.4004e-1, 0.1142e-1, 0.2, 0.4912e1, 0.1428e1), nrow=2, ncol=3, byrow=TRUE) ## Initialize the W matrix
    print('Initial W matrix (before total pivoting):')
    print(W)

    n = dim(W)[1]

    pivot_benefit = abs(W[1:n, 1:n])
    print('pivot_benefit=')
    print(pivot_benefit)

    ## Based on the above we perform elimination of the variable 2 in the non-pivot rows using the row 2:
    ##
    ps_kk = 2
    qs_kk = 2

    ri = 1 
    m = W[ri, qs_kk] = signif(W[ri, qs_kk] / W[ps_kk, qs_kk], n_digits)
    W[ri, 1] = signif(W[ri, 1] - m * W[ps_kk, 1], n_digits) ## <- do the variables for this row
    W[ri, n+1] = signif(W[ri, n+1] - m * W[ps_kk, n+1], n_digits) ## <- do the right-hand-side for this row
    ##print('Final W before backsubstitution:')
    ##print(W)

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

    ## Backsubstitution:
    ##
    kk = 1
    x_c[kk] = W[1, n+1] / W[kk, kk]

    kk = 2
    known_x = x_c[1]
    x_coefficients = W[2, 1]
    x_c[kk] = ( W[kk, n+1] - sum(x_coefficients * known_x) ) / W[kk, kk]

    print(sprintf('Using total pivoting (n_digits in signif= %d); x_c=', n_digits))
    print(x_c)
    print(sprintf('error in x_c= %10.6f', sum((x_c-x_true)^2)))
}