コレスキー分解

Rでコレスキー分解の関数は"chol"
この中身は、と思ったが・・・

Rで
>getS3method("chol","default")
とやって、
RからFortranの関数記載を呼び出し、

            • -

.Fortran("chol", x = x, n, n, v = matrix(0, nrow = n, ncol = n), info = integer(1), DUP = FALSE, PACKAGE = "base")

              • -

という記述があって、

                  • -

>help(.Fortran)

                • -

から、

                  • -

Functions to make calls to compiled code that has been loaded into R.

                  • -

というヘルプ記事があって、Fortranの関数を呼んでいることがわかった
けれども、ここで呼ばれているFortranっていうのは、どこのどれ?



ここにもあるように、そのほかの計算ライブラリを入れておけば、Rから使える
らしいこともわかったものの、それ以上のことがわからなかった。が・・・

追記
RはもともとCとFortranのつぎはぎで作られたソフトなので,Fotranで書かれた関数群がいたるところに点在しているという。
その大部分は R-2.?.?/bin の下の R.dll(おそらく) にアーカイブされていて必要になると関数.Fortranや.CなどによってR側から動的に読み込まれるらしい。

ちなみに関数cholは,Fortranで実装された数値計算ライブラリLINPACKの一つの関数を呼び出している、と、こういう仕組みらしい。

                              • -
              • -

>getS3method("chol","default")
function (x, pivot = FALSE, LINPACK = pivot, ...)
{
if (is.complex(x))
stop("complex matrices not permitted at present")
else if (!is.numeric(x))
stop("non-numeric argument to 'chol'")
if (is.matrix(x)) {
if (nrow(x) != ncol(x))
stop("non-square matrix in 'chol'")
n <- nrow(x)
}
else {
if (length(x) != 1)
stop("non-matrix argument to 'chol'")
n <- as.integer(1)
}
if (!pivot && !LINPACK)
return(.Call("La_chol", as.matrix(x), PACKAGE = "base"))
if (!is.double(x))
storage.mode(x) <- "double"
if (pivot) {
xx <- x
xx[lower.tri(xx)] <- 0
z <- .Fortran("dchdc", x = xx, n, n, double(n), piv = as.integer
(rep.int(0,
n)), as.integer(pivot), rank = integer(1), DUP = FALSE,
PACKAGE = "base")
if (z$rank < n)
if (!pivot)
stop("matrix not positive definite")
else warning("matrix not positive definite")
robj <- z$x
if (pivot) {
attr(robj, "pivot") <- z$piv
attr(robj, "rank") <- z$rank
if (!is.null(cn <- colnames(x)))
colnames(robj) <- cn[z$piv]
}
robj
}
else {
z <- .Fortran("chol", x = x, n, n, v = matrix(0, nrow = n,
ncol = n), info = integer(1), DUP = FALSE, PACKAGE = "base")
if (z$info)
stop("non-positive definite matrix in 'chol'")
z$v
}
}