[F#] 行列式を計算する

F# MathProvider だと、あっさりと MathProvider.LinearAlgebra.inv なる関数があるのだが、自作してみる。

let det ( A' : matrix ) =
    let A = A'.Copy()
    let n = A.NumCols-1
    let B = Matrix.identity A.NumCols
    let mutable det = 1.0

    for k in 0..n do
        // akkを1にする
        let p = A.[k,k]
        det <- det * p
        for j in 0..n do
            A.&#91;k,j&#93; <- A.&#91;k,j&#93; / p
            B.&#91;k,j&#93; <- B.&#91;k,j&#93; / p
        // k行以外から引く
        for i in 0..n do
            if i <> k then
                let d = A.[i,k]
                for j in 0..n do
                    A.[i,j] <- A.&#91;i,j&#93; - A.&#91;k,j&#93; * d
                    B.&#91;i,j&#93; <- B.&#91;i,j&#93; - B.&#91;k,j&#93; * d
    det
&#91;/code&#93;
<p>
逆行列の計算
<a href="http://www.asahi-net.or.jp/~uc3k-ymd/Lesson/Section03/invmat.html">http://www.asahi-net.or.jp/~uc3k-ymd/Lesson/Section03/invmat.html</a>
</p>
<p>
の中にある Fortran のコードから擬似写しなんだけど、「det <- det * p」なところで、行列式を保存。他はいらないはずなんだが、うまく括りだせないのでそのままです。
数学的には det A = 0 の場合は、逆行列が計算できないので連立一次方程式が解けないことなるのだが、有限要素法の場合には拘束条件を間違えない限り、逆行列は存在するので問題なし（なはず）。この方式だと det の計算自体に逆行列の計算が含まれてしまっているので本末転倒だし、計算時間もかかってしまう。なので、実際には 0 除算にならなければ逆行列の計算を続ける、という方式でよいのかも。
</p>
[code lang="csharp"]
let A = matrix [[3.0; 1.0;]; 
                [2.0; 5.0;]]

det A

な感じで計算ができて、

val it : float = 13.0

の結果が得られる。