linear algebra - How to do Gaussian elimination in R (do not use "solve")

Question

A <- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
U.pls <- cbind(A,b)
for (i in 1:p){
    for (j in (i+1):(p+1)) U.pls[i,j] <- U.pls[i,j]/U.pls[i,i]
    U.pls[i,i] <- 1
    if (i < p) {
       for (k in (i+1):p) U.pls[k,] <- 
                   U.pls[k,] - U.pls[k,i]/U.pls[i,i]*U.pls[i,]
    }
}
U.pls
x <- rep(0,p)
for (i in p:1){
  if (i < p){
     temp <- 0
     for (j in (i+1):p) temp <- temp + U.pls[i,j]*x[j]
     x[i] <- U.pls[i,p+1] - temp
  }
  else x[i] <- U.pls[i,p+1]
}
x

output

> U.pls
     [,1] [,2] [,3] [,4]
[1,]    1 -2.5    2 -1.5
[2,]    0  1.0 -Inf  Inf
[3,]    0  0.0    1  NaN
> x
[1] NaN NaN NaN

In this way, I cannot solve solution in some case. Even though I know the reason why the error happens mathematically, I cannot fix the error in R. Help me with some fix. Thank you in advance.

Answer 1

I had to reorder you matrix A, don't know how do a generic code:

A <- matrix(c(2,-5,4,1,-4,6,1,-2.5,1),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

for (i in 2:p){
 for (j in i:p) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
}

for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls

EDIT:

A <- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
b <- matrix(c(-3,5,10),nrow=3,ncol=1)
p <- nrow(A)
(U.pls <- cbind(A,b))

U.pls[1,] <- U.pls[1,]/U.pls[1,1]

i <- 2
while (i < p+1) {

 j <- i
 while (j < p+1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i-1, ] * U.pls[j, i-1]
  j <- j+1
 }
 while (U.pls[i,i] == 0) {
  U.pls <- rbind(U.pls[-i,],U.pls[i,])
 }
 U.pls[i,] <- U.pls[i,]/U.pls[i,i]
 i <- i+1
}
for (i in p:2){
 for (j in i:2-1) {
  U.pls[j, ] <- U.pls[j, ] - U.pls[i, ] * U.pls[j, i]
 }
}
U.pls