Мне пришлось изменить порядок вашей матрицы A, не знаю, как сделать общий код:

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

РЕДАКТИРОВАТЬ:

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

Извините, если этот ответ запоздал.
Я написал функцию для исключения Гаусса.

Примечание:
i) Эта функция работает только с действительными числами, а не с переменными.
ii) В ответе в конце я использовал функцию as.fractions(). Она доступна только в пакете MASS и вам нужно иметь как минимум R версии 3.2.5 или новее, чтобы использовать его. Вы можете опустить это, если хотите, но элементы матриц будут десятичными.
iii) Эта функция также возвращает матрицы, используемые в разложении PA = LDU Если вы хотите использовать только разложение PA = LU в качестве второго аргумент ЛОЖЬ.

Надеюсь, это поможет!

gauss_elimination = function(matrix,diagonal = T){
  #This function performs the gaussian elimination for one column
  for_one_column = function(matrix, starting_row = 1, column = 1){
    for (i in (starting_row + 1):nrow(matrix)){
      L_table[i,column] <<-  matrix[i,column]/matrix[starting_row,column]
      matrix[i,] = matrix[i,] - matrix[i,column]/
        matrix[starting_row,column]*
        matrix[starting_row,]
    }
    return(matrix)
  }

  #This function changes a row of a matrix with another one
  change_lines = function(matrix,line_1,line_2,column_1 = 1,
                          column_2 =  ncol(matrix)){
    scapegoat = matrix[line_1,column_1:column_2]
    matrix[line_1,column_1:column_2] = matrix[line_2,column_1:column_2]
    matrix[line_2,column_1:column_2] = scapegoat
    return(matrix)
  }

  #This function checks if alternations need to be made
  checking_zeroes = function(matrix,starting_row,column){
    #If pilot is not zero
    if (matrix[starting_row,column] != 0){
      return ("No alternation needed")
      #If pilot is zero
    }else{
      row_to_change = 0
      for (i in starting_row:nrow(matrix)){
        if (matrix[i,column] != 0){
          row_to_change = i
          break
        }
      }
      #If both the pilot and all the elements below it are zero
      if (row_to_change == 0){
        return("Skip this column")
        #If the pilot is zero but at least one below it is not
      }else{
        return(c("Alternation needed",row_to_change))
      }
    }
  }

  #The main program
  row_to_work = 1 ; alternation_table = diag(nrow(matrix))
  L_table = diag(nrow(matrix))
  for (column_to_work in 1:ncol(matrix)){
    a = checking_zeroes(matrix,row_to_work,column_to_work)
    if (a[1] == "Alternation needed"){
      matrix = change_lines(matrix,as.numeric(a[2]),row_to_work)
      alternation_table = change_lines(alternation_table,row_to_work,
                                       as.numeric(a[2]))
      if (as.numeric(a[2]) != 1){
        L_table = change_lines(L_table,row_to_work,
                               as.numeric(a[2]),1,column_to_work - 1)
      }
    }
    if (a[1] == "Skip this column"){
      next()
    }
    matrix = for_one_column(matrix,row_to_work,column_to_work)
    if(row_to_work + 1 == nrow(matrix)){
      break
    }
    row_to_work = row_to_work + 1
  }
  if (diagonal == FALSE){
    return(list("P" = as.fractions(alternation_table),
                "L" = as.fractions(L_table),
                "U : The matrix after the elimination"=as.fractions(matrix)))
  }else{
    D = diag(nrow(matrix))
    diag(D) = diag(matrix)
    for (i in 1:nrow(D)){
      matrix[i,] = matrix[i,]/diag(D)[i]
    }
    return(list("P" = as.fractions(alternation_table),
                "L" = as.fractions(L_table),
                "D" = as.fractions(D),
                "U : The matrix after the elimination"=as.fractions(matrix)))
  }
}

‹-функция Гаусса () { A ‹- матрица ( c ( 2 , -5 , 4 , 1 , 4 , 1 , 1 , -4 , 6 ), byrow = T , nrow = 3 , ncol = 3 ) b ‹- матрица ( c ( -3 , 5 , 10 ), nrow = 3 , ncol = 1 ) U.pls ‹- cbind (A,b) p ‹- nrow ( A )

  r <- ncol ( U.pls )
  X <- matrix ( c(rep (0,p)))
  for ( i in 1 :(p-1))  {
    cons <-    U.pls[i+1,i]/ U.pls[i,i]
    for ( j in (i+1):(p))   {
        U.pls[j,] <-  U.pls[j,] -   U.pls[i,]  *   cons
            }
     }

    X [p] <- U.pls[ p,r]/ U.pls[p,p]

  for (k in (p-1):1)      {
       suma = 0
       for (l in (k+1):p) {
            suma = suma + U.pls[k,l]*  X[l]  }
          X[k] <- (U.pls[k,r]-   suma   )/ U.pls[k,k]
          }
  return (X)

}

Вы можете выполнить исключение Гаусса через библиотеку pracma.

Чтобы установить его, вы можете ввести:

install.packages("pracma")

Затем используйте его следующим образом:

library(pracma)

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)

rref(cbind(A, b)

Результат:

     [,1] [,2] [,3]  [,4]
[1,]    1    0    0 -51.0
[2,]    0    1    0 -25.0
[3,]    0    0    1  -6.5