R で書いたガウスの消去法

#
#	ガウスの消去法
# 	...modified on 2015.05.18 (04:02)
#
# この中身を I:￥rdata￥gauss.r に保存した上で、
# source('I:/rdata/gauss.r')
# m <- gauss(matrix(c(1,2,3,2,3,4,3,4,6,2,4,5),3,4)) 
# cat(m) のようにして使います
#
N <- 3
abs <- function(f) {ifelse(f<0,-f,f)}
gauss <- function(a) {
  # 前処理: 行を大きい順に並び替える
  for (i in 1:N) {
    max <- 0
    pvt <- i
    for (j in i:N) {
      if (abs(a[j,i]) > max) {
        max <- abs(a[j,i])
        pvt <- j
      }
    }
    if (pvt != i) {
      for (j in 1:(N+1)) {
      	tmp <- a[i,j]
	a[i,j] <- a[pvt,j]
	a[pvt,j] <- tmp
      }
    }
  }
  # 本処理
  for (i in 1:N) {
    p <- a[i,i] # 対角要素を保存して
    a[i,i] <- 1 # 対角要素を１にする
    for (j in (i+1):(N+1)) {
      a[i,j] <- a[i,j]/p
    }
    for (j in (i+1):N) {
      if (j>N) break
      q <- a[j,i]
      for (k in (i+1):(N+1)) {
        a[j,k] <- a[j,k] - q * a[i,k]
      }
      a[j,i] <- 0 # zero 要素
    }
  }
  # 解を計算する
  x <- array(0,N)
  for (i in N:1) {
     x[i] = a[i,N+1]
     for (j in N:(i+1)) {
       if ((i+1)>N) next
       x[i] <- x[i] - a[i,j]*x[j]
     }
  }
  return(x)
}