LU分解を F# らしく書き直してみます。
せっかくの F# なのに for ループを使っているのがもったいないし（？）、数式のΣから外れています。ここを sum 関数を使って書き直したいのですが… matrix に対しての Seq.sum とかは要素全体に掛け算してしまうので無理（だと思う）なので、
C#er のためのやさしい再帰入門: いげ太のブログ
http://igeta.cocolog-nifty.com/blog/2011/02/tailcall.html#more
を参考にして、再帰関数化しています。
/// LU分解
let LU2 ( A : matrix ) =
let n = A.NumCols
let L = Matrix.identity n
let U = Matrix.zero n n
// sum(k=0,n) L(i,k)*U(k,j) の計算を関数化
// 余計わかりづらいか？
let sum i j k =
let rec f i j k acc =
if k < 0 then 0.0
elif k > 0 then f i j (k-1) (acc + L.[i,k]*U.[k,j])
else L.[i,k]*U.[k,j]
f i j k 0.0
for j in 0..n-1 do
// Uを計算
for i in 0..j do
U.[i,j] <- A.[i,j] - sum i j (i-1)
// Lを計算
for i in j+1..n-1 do
L.[i,j] <- (A.[i,j] - sum i j (j-1))/U.[j,j]
// 結果を返す
(L,U)
[/code]
LとUを計算するところで、ΣL.[i,k]*U.[k.j] が同じなので、括りだしてみたのですが、余計わかりづらいかも。最後の for ループに適用するとよいのかもしれません。まあ、この程度の簡単なものならば再帰を使うとかえってややこしいという感じですかね。
ならば、for の二重ループの中で、sum を定義して、ちょっと簡単にしてみたのが次のコードです。
/// LU分解
let LU3 ( A : matrix ) =
let n = A.NumCols
let L = Matrix.identity n
let U = Matrix.zero n n
// sum(k=0,n) L(i,k)*U(k,j) の計算を関数化
// 普通にforループで合計を計算する
for j in 0..n-1 do
// Uを計算
for i in 0..n-1 do
let sum max =
let mutable s = 0.0
for k in 0..max do
s <- s + L.[i,k]*U.[k,j]
s
if i <= j then
// Uを計算
U.[i,j] <- A.[i,j] - sum (i-1)
else
// Lを計算
L.[i,j] <- (A.[i,j] - sum (j-1))/U.[j,j]
// 結果を返す
(L,U)
[/code]
ループの中で関数が定義できるから、便利ってのがありますね。ループ変数のi,jもそのまま使えるので、sum 関数の引数が少なくて済みます。
お次は、このLU分解を使って、二次元トラスを具体的に解いてみます。
