コレスキー法

home > さまざまな言語で数値計算 > 連立一次方程式 >

コレスキー法

コレスキー法を使って, 連立一次方程式の解を求める .

行と列とを入れ替えても(転置行列)一致する行列を対称行列と言う.

コレスキー法では, 対称行列 $A$ を下三角行列 $L$ と, $L$ の転置行列 $L^T$ との積に分解し $ ( A = LL^T ) $ ,

$ Ly = b $ から $y$ を求め, $ L^Tx = y $ から $x$ を求める.

VBScript

Option Explicit

Private Const N = 4

Private a: a = Array(Array(5.0,2.0,3.0,4.0),Array(2.0,10.0,6.0,7.0),Array(3.0,6.0,15.0,9.0),Array(4.0,7.0,9.0,20.0))
Private b: b = Array(34.0,68.0,96.0,125.0)

WScript.StdOut.WriteLine "A"
disp_matrix a
WScript.StdOut.WriteLine "B"
disp_vector b
WScript.StdOut.WriteLine

'前進消去
forward_elimination a,b

WScript.StdOut.WriteLine "LL^T"
disp_matrix a

'Ly=b から y を求める (前進代入)
Private y: y = Array(0.0,0.0,0.0,0.0)
forward_substitution a,b,y

WScript.StdOut.WriteLine "Y"
disp_vector y

'L^Tx=y から x を求める (後退代入)
Private x: x = Array(0.0,0.0,0.0,0.0)
backward_substitution a,y,x

WScript.StdOut.WriteLine "X"
disp_vector x

'前進消去
Private Sub forward_elimination(ByRef a, ByVal b)
    Dim pivot, row, col
    For pivot = 0 To (N - 1)
        Dim s: s = 0
        For col = 0 To (pivot - 1)
            s = s + a(pivot)(col) * a(pivot)(col)
        Next
        'ここで根号の中が負の値になると計算できない！
        a(pivot)(pivot) = Sqr(a(pivot)(pivot) - s)

        For row = (pivot + 1) To (N - 1)
            s = 0
            For col = 0 To (pivot - 1)
                s = s + a(row)(col) * a(pivot)(col)
            Next
            a(row)(pivot) = (a(row)(pivot) - s) / a(pivot)(pivot)
            a(pivot)(row) = a(row)(pivot)
        Next
    Next
End Sub
'前進代入
Private Sub forward_substitution(ByVal a, ByVal b, ByRef y)
    Dim row, col
    For row = 0 To (N - 1)
        For col = 0 To row
            b(row) = b(row) - a(row)(col) * y(col)
        Next
        y(row) = b(row) / a(row)(row)
    Next
End Sub
'後退代入
Private Sub backward_substitution(ByVal a, ByVal y, ByRef x)
    Dim row, col
    For row = (N - 1) To 0 Step -1
        For col = (N - 1) To (row + 1) Step - 1
            y(row) = y(row) - a(row)(col) * x(col)
        Next
        x(row) = y(row) / a(row)(row)
    Next
End Sub
'１次元配列を表示
Private Sub disp_vector(ByVal row())
    Dim col
    For Each col In row
        WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab
    Next
    WScript.StdOut.WriteLine
End Sub
'２次元配列を表示
Private Sub disp_matrix(ByVal matrix)
    Dim row, col
    For Each row In matrix
        For Each col In row
            WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab
        Next
        WScript.StdOut.WriteLine
    Next
End Sub

Z:\>cscript //nologo Z:\1006.vbs
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

JScript

var N = 4

var a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]]
var b = [34,68,96,125]

WScript.StdOut.WriteLine("A")
disp_matrix(a)
WScript.StdOut.WriteLine("B")
disp_vector(b)
WScript.StdOut.WriteLine()

// 前進消去
forward_elimination(a,b)

WScript.StdOut.WriteLine("LL^T")
disp_matrix(a)

// Ly=b から y を求める (前進代入)
var y = [0,0,0,0]
forward_substitution(a,b,y)

WScript.StdOut.WriteLine("Y")
disp_vector(y)

// L^Tx=y から x を求める (後退代入)
var x = [0,0,0,0]
backward_substitution(a,y,x)

WScript.StdOut.WriteLine("X")
disp_vector(x)

// 前進消去
function forward_elimination(a, b)
{
    for (pivot = 0; pivot < N; pivot++)
    {
        var s = 0
        for (col = 0; col < pivot; col++)
            s += a[pivot][col] * a[pivot][col]
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s)

        for (row = pivot + 1; row < N; row++)
        {
            s = 0
            for (col = 0; col < pivot; col++)
                s += a[row][col] * a[pivot][col]
            a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot]
            a[pivot][row] =  a[row][pivot]
        }
    }
}
// 前進代入
function forward_substitution(a, b, y)
{
    for (row = 0; row < N; row++)
    {
        for (col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col]
        y[row] = b[row] / a[row][row]
    }
}
// 後退代入
function backward_substitution(a, y, x)
{
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col]
        x[row] = y[row] / a[row][row]
    }
}
// １次元配列を表示
function disp_vector(row)
{
    for (var i = 0; i < N; i++)
        WScript.StdOut.Write(("              "    + row[i].toFixed(10)).slice(-14) + "\t")
    WScript.StdOut.WriteLine()
}
// ２次元配列を表示
function disp_matrix(matrix)
{
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
            WScript.StdOut.Write(("              "    + matrix[i][j].toFixed(10)).slice(-14) + "\t")
        WScript.StdOut.WriteLine()
    }
}

Z:\>cscript //nologo Z:\1006.js
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

PowerShell

set-variable -name N -value 3 -option constant

# 前進消去
function forward_elimination($a, $b)
{
    foreach ($pivot in 0..$N)
    {
        $s = 0
        for ($col = 0; $col -lt $pivot; $col++)
        {
            $s += $a[$pivot][$col] * $a[$pivot][$col]
        }
        # ここで根号の中が負の値になると計算できない！
        $a[$pivot][$pivot] = [Math]::Sqrt($a[$pivot][$pivot] - $s)

        for ($row = $pivot + 1; $row -le $N; $row++)
        {
            $s = 0
            for ($col = 0; $col -lt $pivot; $col++)
            {
                $s += $a[$row][$col] * $a[$pivot][$col]
            }
            $a[$row][$pivot] = ($a[$row][$pivot] - $s) / $a[$pivot][$pivot]
            $a[$pivot][$row] = $a[$row][$pivot]
        }
    }
}
# 前進代入
function forward_substitution($a, $b, $y)
{
    foreach ($row in 0..$N)
    {
        foreach ($col in 0..($row - 1))
        {
            $b[$row] -= $a[$row][$col] * $y[$col]
        }
        $y[$row] = $b[$row] / $a[$row][$row]
    }
}
# 後退代入
function backward_substitution($a, $y, $x)
{
    foreach ($row in $N..0)
    {
        if (($row+1) -le $N)
        {
            foreach ($col in $N..($row + 1))
            {
                $y[$row] -= $a[$row][$col] * $x[$col]
            }
        }
        $x[$row] = $y[$row] / $a[$row][$row]
    }
}
# １次元配列を表示
function disp_vector($row)
{
    foreach ($col in $row)
    {
        Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline
    }
    Write-Host
}
# ２次元配列を表示
function disp_matrix($matrix)
{
    foreach ($row in $matrix)
    {
        foreach ($col in $row)
        {
            Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline
        }
        Write-Host
    }
}

$a = (5.0, 2.0, 3.0, 4.0), (2.0, 10.0, 6.0, 7.0), (3.0, 6.0, 15.0, 9.0), (4.0, 7.0, 9.0, 20.0)
$b = 34.0, 68.0, 96.0, 125.0

Write-Host "A"
disp_matrix $a
Write-Host "B"
disp_vector $b
Write-Host ""

# 前進消去
forward_elimination $a $b

Write-Host "LL^T"
disp_matrix $a

# Ly=b から y を求める (前進代入)
$y = 0.0, 0.0, 0.0, 0.0
forward_substitution $a $b $y

Write-Host "Y"
disp_vector $y

# L^Tx=y から x を求める (後退代入)
$x = 0.0, 0.0, 0.0, 0.0
backward_substitution $a $y $x

Write-Host "X"
disp_vector $x

Z:\>powershell -file Z:\1006.ps1
A
  5.0000000000  2.0000000000  3.0000000000  4.0000000000
  2.0000000000 10.0000000000  6.0000000000  7.0000000000
  3.0000000000  6.0000000000 15.0000000000  9.0000000000
  4.0000000000  7.0000000000  9.0000000000 20.0000000000
B
 34.0000000000 68.0000000000 96.0000000000125.0000000000

LL^T
  2.2360679775  0.8944271910  1.3416407865  1.7888543820
  0.8944271910  3.0331501776  1.5825131361  1.7803272782
  1.3416407865  1.5825131361  3.2704207946  1.1566122322
  1.7888543820  1.7803272782  1.1566122322  3.5060922587
Y
 15.2052622470 17.9351488764 14.4377113129 14.0243690350
X
  1.0000000000  2.0000000000  3.0000000000  4.0000000000

Perl

use constant N => 3;

my @a = ([5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]);
my @b = (34,68,96,125);

print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 前進消去
forward_elimination(\@a,\@b);

print "LL^T\n";
disp_matrix(\@a);

# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);

print "Y\n";
disp_vector(\@y);

# L^Tx=y から x を求める (後退代入)
my @x = (0,0,0,0);
backward_substitution(\@a,\@y,\@x);

print "X\n";
disp_vector(\@x);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        my $s = 0;
        for $col (0..($pivot-1))
        {
            $s += $$a[$pivot][$col] * $$a[$pivot][$col];
        }
        # ここで根号の中が負の値になると計算できない！
        $$a[$pivot][$pivot] = sqrt($$a[$pivot][$pivot] - $s);

        for $row (($pivot + 1)..N)
        {
            $s = 0;
            for $col (0..($pivot-1))
            {
                $s += $$a[$row][$col] * $$a[$pivot][$col];
            }
            $$a[$row][$pivot] = ($$a[$row][$pivot] - $s) / $$a[$pivot][$pivot];
            $$a[$pivot][$row] = $$a[$row][$pivot];
        }
    }
}
# 前進代入
sub forward_substitution
{
    my ($a, $b, $y) = @_;

    for $row (0..N)
    {
        for $col (0..$row)
        {
            $$b[$row] -= $$a[$row][$col] * $$y[$col];
        }
        $$y[$row] = $$b[$row] / $$a[$row][$row];
    }
}
# 後退代入
sub backward_substitution
{
    my ($a, $y, $x) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$y[$row] -= $$a[$row][$col] * $$x[$col];
        }
        $$x[$row] = $$y[$row] / $$a[$row][$row];
    }
}
# ２次元配列を表示
sub disp_matrix
{
    my ($matrix) = @_;
    foreach $row (@$matrix)
    {
        foreach $col (@$row)
        {
            printf("%14.10f\t", $col);
        }
        print "\n";
    }
}
# １次元配列を表示
sub disp_vector
{
    my ($row) = @_;
    foreach $col (@$row)
    {
        printf("%14.10f\t", $col);
    }
    print "\n";
}

Z:\>perl Z:\1006.pl
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

PHP

<?php
define("N", 3);

$a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]];
$b = [34,68,96,125];

print "A\n";
disp_matrix($a);
print "B\n";
disp_vector($b);
print "\n";

# 前進消去
forward_elimination($a,$b);

print "LL^T\n";
disp_matrix($a);

# Ly=b から y を求める (前進代入)
$y = [0,0,0,0];
forward_substitution($a,$b,$y);

print "Y\n";
disp_vector($y);

# L^Tx=y から x を求める (後退代入)
$x = [0,0,0,0];
backward_substitution($a,$y,$x);

print "X\n";
disp_vector($x);

# 前進消去
function forward_elimination(&$a, $b)
{
    foreach (range(0, N) as $pivot)
    {
        $s = 0;
        for ($col = 0; $col < $pivot; $col++)
            $s += $a[$pivot][$col] * $a[$pivot][$col];
        # ここで根号の中が負の値になると計算できない！
        $a[$pivot][$pivot] = sqrt($a[$pivot][$pivot] - $s);

        for ($row = $pivot + 1; $row <= N; $row++)
        {
            $s = 0;
            for ($col = 0; $col < $pivot; $col++)
                $s += $a[$row][$col] * $a[$pivot][$col];
            $a[$row][$pivot] = ($a[$row][$pivot] - $s) / $a[$pivot][$pivot];
            $a[$pivot][$row] = $a[$row][$pivot];
        }
    }
}
# 前進代入
function forward_substitution($a, $b, &$y)
{
    foreach (range(0, N) as $row)
    {
        if ($row != 0)
            foreach (range(0, $row - 1) as $col)
                $b[$row] -= $a[$row][$col] * $y[$col];
        $y[$row] = $b[$row] / $a[$row][$row];
    }
}
# 後退代入
function backward_substitution($a, $y, &$x)
{
    foreach (range(N, 0) as $row)
    {
        if (($row + 1) <= N)
            foreach (range(N, $row + 1) as $col)
                $y[$row] -= $a[$row][$col] * $x[$col];
        $x[$row] = $y[$row] / $a[$row][$row];
    }
}
# ２次元配列を表示
function disp_matrix($matrix)
{
    foreach ($matrix as $row)
    {
        foreach ($row as $col)
            printf("%14.10f\t", $col);
        print "\n";
    }
}
# １次元配列を表示
function disp_vector($row)
{
    foreach ($row as $col)
        printf("%14.10f\t", $col);
    print "\n";
}
?>

Z:\>php Z:\1006.php
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Python

# coding: Shift_JIS

import math

N = 4

# 前進消去
def forward_elimination(a, b):
    for pivot in range(0, N, 1):
        s = 0.0
        for col in range(0, pivot, 1):
            s += a[pivot][col] * a[pivot][col]
        a[pivot][pivot] = math.sqrt(a[pivot][pivot] - s)

        for row in range(pivot + 1, N, 1):
            s = 0.0
            for col in range(0, pivot, 1):
                s += a[row][col] * a[pivot][col]
            a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot]
            a[pivot][row] = a[row][pivot]

# 前進代入
def forward_substitution(a, b, y):
    for row in range(0, N, 1):
        for col in range(0, row + 1, 1):
            b[row] -= a[row][col] * y[col]

        y[row] = b[row] / a[row][row]

# 後退代入
def backward_substitution(a, y, x):
    for row in range(N - 1, -1, -1):
        for col in range(N - 1, row, -1):
            y[row] -= a[row][col] * x[col]
        x[row] = y[row] / a[row][row]

# ２次元配列を表示
def disp_matrix(matrix):
    for row in matrix:
        for col in row:
            print "%14.10f\t" % col,
        print ""

# １次元配列を表示
def disp_vector(row):
    for col in row:
        print "%14.10f\t" % col,
    print ""

a = [[5.0,2.0,3.0,4.0],[2.0,10.0,6.0,7.0],[3.0,6.0,15.0,9.0],[4.0,7.0,9.0,20.0]]
b =  [34.0,68.0,96.0,125.0]

print "A"
disp_matrix(a)
print "B"
disp_vector(b)
print ""

# 前進消去
forward_elimination(a,b)

print "LL^T"
disp_matrix(a)

# Ly=b から y を求める (前進代入)
y = [0.0, 0.0, 0.0, 0.0]
forward_substitution(a,b,y)

print "Y"
disp_vector(y)

# L^Tx=y から x を求める (後退代入)
x =  [ 0.0,  0.0,  0.0,  0.0]
backward_substitution(a,y,x)

print "X"
disp_vector(x)

Z:\>python Z:\1006.py
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Ruby

# coding: Shift_JIS

N = 3

# ２次元配列を表示
def disp_matrix(matrix)
    matrix.each do |row|
        row.each do |col|
            printf("%14.10f\t", col)
        end
        puts ""
    end
end
# １次元配列を表示
def disp_vector(row)
    row.each do |col|
        printf("%14.10f\t", col)
    end
    puts ""
end
# 前進消去
def forward_elimination(a, b)
    0.upto(N) do |pivot|
        s = 0
        0.upto(pivot-1) do |col|
            s += a[pivot][col] * a[pivot][col]
        end
        # ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s)

        (pivot + 1).upto(N) do |row|
            s = 0
            0.upto(pivot-1) do |col|
                s += a[row][col] * a[pivot][col]
            end
            a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot]
            a[pivot][row] = a[row][pivot]
        end
    end
end
# 前進代入
def forward_substitution(a, b, y)
    (0..N).each do |row|
        (0..row).each do |col|
            b[row] -= a[row][col] * y[col]
        end
        y[row] = b[row] / a[row][row]
    end
end
# 後退代入
def backward_substitution(a, y, x)
    N.downto(0) do |row|
        N.downto(row + 1) do |col|
            y[row] -= a[row][col] * x[col]
        end
        x[row] = y[row] / a[row][row]
    end
end

a = [[5.0,2.0,3.0,4.0],[2.0,10.0,6.0,7.0],[3.0,6.0,15.0,9.0],[4.0,7.0,9.0,20.0]]
b =  [34.0,68.0,96.0,125.0]

puts "A"
disp_matrix(a)
puts "B"
disp_vector(b)
puts ""

# 前進消去
forward_elimination(a,b)

puts "LL^T"
disp_matrix(a)

# Ly=b から y を求める (前進代入)
y = [0.0,0.0,0.0,0.0]
forward_substitution(a,b,y)

puts "Y"
disp_vector(y)

# L^Tx=y から x を求める (後退代入)
x =  [0.0,0.0,0.0,0.0]
backward_substitution(a,y,x)

puts "X"
disp_vector(x)

Z:\>ruby Z:\1006.rb
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Groovy

N = 3

def a = [[5.0,2.0,3.0,4.0],[2.0,10.0,6.0,7.0],[3.0,6.0,15.0,9.0],[4.0,7.0,9.0,20.0]]  as double[][]
def b =  [34.0,68.0,96.0,125.0] as double[]

println("A")
disp_matrix(a)
println("B")
disp_vector(b)
println()

// 前進消去
forward_elimination(a,b)

println("LL^T")
disp_matrix(a)
println()

// Ly=b から y を求める (前進代入)
y = [0.0,0.0,0.0,0.0]
forward_substitution(a,b,y)

println "Y"
disp_vector(y)

// L^Tx=y から x を求める (後退代入)
x =  [0.0,0.0,0.0,0.0]
backward_substitution(a,y,x)

println("X")
disp_vector(x)

// 前進消去
def forward_elimination(a, b) {
    for (pivot in 0..N) {
        s = 0
        if (pivot > 0)
            for (col in 0..(pivot-1))
                s += a[pivot][col] * a[pivot][col]
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s)

        if (pivot >= N) continue

        for (row in (pivot+1)..N) {
            s = 0
            if (pivot > 0)
                for (col in 0..(pivot-1))
                    s += a[row][col] * a[pivot][col]
            a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot]
            a[pivot][row] = a[row][pivot]
        }
    }
}
// 前進代入
def forward_substitution(a, b, y) {
    for (row in 0..N) {
        for (col in 0..row)
            b[row] -= a[row][col] * y[col]
        y[row] = b[row] / a[row][row]
    }
}
// 後退代入
def backward_substitution(a, y, x) {
    for (row = N; row >= 0; row--) {
        for (col = N; col > row; col--)
            y[row] -= a[row][col] * x[col]
        x[row] = y[row] / a[row][row]
    }
}
// １次元配列を表示
def disp_matrix(matrix) {
    for (row in matrix) {
        for (col in row) {
            printf ("%14.10f\t" , col)
        }
        println()
    }
}
// ２次元配列を表示
def disp_vector(row) {
    for (col in row) {
        printf ("%14.10f\t" , col)
    }
    println()
}

Z:\>groovy Groovy1006.groovy
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587

Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Pascal

program Pas1006(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    N = 3;

type
    TwoDimArray = array [0..N,0..N] of Double;

// １次元配列を表示
procedure disp_vector(row:array of Double);
var
    i:Integer;
begin
    for i := Low(row) to High(row) do
        write(format('%14.10f'#9, [row[i]]));
    writeln();
end;
// ２次元配列を表示
procedure disp_matrix(matrix:TwoDimArray);
var
    i,j:Integer;
begin
    for i := Low(matrix) to High(matrix) do
    begin
        for j := Low(matrix) to High(matrix) do
            write(format('%14.10f'#9, [matrix[i,j]]));
        writeln();
    end;
end;
// 前進消去
procedure forward_elimination(var a:TwoDimArray; var b:array of Double);
var
    s:Double;
    pivot, row, col:Integer;
begin
    for pivot := Low(a) to High(a) do
    begin
        s := 0;
        for col := Low(a) to pivot - 1 do
            s := s + a[pivot,col] * a[pivot,col];
        // ここで根号の中が負の値になると計算できない！
        a[pivot,pivot] := Sqrt(a[pivot,pivot] - s);

        for row := (pivot + 1) to High(a) do
        begin
            s := 0;
            for col := Low(a) to pivot - 1 do
                s := s + a[row,col] * a[pivot,col];
            a[row,pivot] := (a[row,pivot] - s) / a[pivot,pivot];
            a[pivot,row] := a[row,pivot];
        end;
    end;
end;
// 前進代入
procedure forward_substitution(a:TwoDimArray; b:array of Double; var y:array of Double);
var
    row, col:Integer;
begin
    for row := Low(a) to High(a) do
    begin
        for col := Low(a) to row do
            b[row] := b[row] - a[row,col] * y[col];
        y[row] := b[row] / a[row,row];
    end;
end;
// 後退代入
procedure backward_substitution(a:TwoDimArray; y:array of Double; var x:array of Double);
var
    row, col:Integer;
begin
    for row := High(a) downto Low(a) do
    begin
        for col := High(a) downto row + 1 do
            y[row] := y[row] - a[row,col] * x[col];
        x[row] := y[row] / a[row,row];
    end;
end;

var
    a :TwoDimArray = ((5.0,2.0,3.0,4.0),(2.0,10.0,6.0,7.0),(3.0,6.0,15.0,9.0),(4.0,7.0,9.0,20.0));
    b :array [0..N] of Double = (34.0,68.0,96.0,125.0);
    x :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0);
    y :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0);
begin
    writeln('A');
    disp_matrix(a);
    writeln('B');
    disp_vector(b);
    writeln();

    // 前進消去
    forward_elimination(a,b);

    writeln('LL^T');
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    forward_substitution(a,b,y);

    writeln('Y');
    disp_vector(y);

    // L^Tx=y から x を求める (後退代入)
    backward_substitution(a,y,x);

    writeln('X');
    disp_vector(x);
end.

Z:\>fpc -v0 -l- Pas1006.pp

Z:\>Pas1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Ada

with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions;
use  TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions;

procedure Ada1006 is

    N : Constant Integer := 3;

    type Long_Float_Array       is array (0..N)       of Long_Float;
    type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float;

    a:Long_Float_TwoDimArray := ((5.0,2.0,3.0,4.0),(2.0,10.0,6.0,7.0),(3.0,6.0,15.0,9.0),(4.0,7.0,9.0,20.0));
    b:Long_Float_Array := (34.0,68.0,96.0,125.0);
    y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0);
    x:Long_Float_Array := (0.0, 0.0, 0.0, 0.0);

    -- １次元配列を表示
    procedure disp_vector(row:Long_Float_Array) is
    begin
        for i in row'Range loop
            Put(row(i), Fore=>3, Aft=>10, Exp=>0);
            Put(Ascii.HT);
        end loop;
        New_Line;
    end disp_vector;

    -- ２次元配列を表示
    procedure disp_matrix(matrix:Long_Float_TwoDimArray) is
    begin
        for i in matrix'Range loop
            for j in matrix'Range loop
                Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0);
                Put(Ascii.HT);
            end loop;
            New_Line;
        end loop;
    end disp_matrix;

    -- 前進消去
    procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is
        s:Long_Float;
    begin
        for pivot in a'Range loop
            s := 0.0;
            for col in a'First..(pivot - 1) loop
                s := s + a(pivot,col) * a(pivot,col);
            end loop;
            -- ここで根号の中が負の値になると計算できない！
            a(pivot,pivot) := Sqrt(a(pivot,pivot) - s);

            for row in (pivot + 1)..a'Last loop
                s := 0.0;
                for col in a'First..(pivot - 1) loop
                    s := s + a(row,col) * a(pivot,col);
                end loop;
                a(row,pivot) := (a(row,pivot) - s) / a(pivot,pivot);
                a(pivot,row) := a(row,pivot);
            end loop;
        end loop;
    end forward_elimination;

    -- 前進代入
    procedure forward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array; y:in out Long_Float_Array) is
    begin
        for row in a'Range loop
            for col in a'First..row loop
                b(row) := b(row) - a(row,col) * y(col);
            end loop;
            y(row) := b(row) / a(row,row);
        end loop;
    end forward_substitution;

    -- 後退代入
    procedure backward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; x:in out Long_Float_Array) is
    begin
        for row in reverse a'Range loop
            for col in reverse (row + 1)..a'Last loop
                y(row) := y(row) - a(row,col) * x(col);
            end loop;
            x(row) := y(row) / a(row,row);
        end loop;
    end backward_substitution;

begin
    Put_Line("A");
    disp_matrix(a);
    Put_Line("B");
    disp_vector(b);
    Put_Line("");

    -- 前進消去
    forward_elimination(a,b);

    Put_Line("LL^T");
    disp_matrix(a);

    -- Ly=b から y を求める (前進代入)
    forward_substitution(a,b,y);

    Put_Line("Y");
    disp_vector(y);

    -- L^Tx=y から x を求める (後退代入)
    backward_substitution(a,y,x);

    Put_Line("X");
    disp_vector(x);
end Ada1006;

xxxxxx@yyyyyy /Z
$ gnatmake Ada1006.adb

xxxxxx@yyyyyy /Z
$ Ada1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

VB.NET

Option Explicit

Module VB1006
    Private Const N As Integer = 3

    Public Sub Main()
        Dim a(,) As Double = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}}
        Dim b()  As Double = {34,68,96,125}

        Console.WriteLine("A")
        disp_matrix(a)
        Console.WriteLine("B")
        disp_vector(b)
        Console.WriteLine

        '前進消去
        forward_elimination(a,b)

        Console.WriteLine("LL^T")
        disp_matrix(a)

        'Ly=b から y を求める (前進代入)
        Dim y()  As Double = {0,0,0,0}
        forward_substitution(a,b,y)

        Console.WriteLine("Y")
        disp_vector(y)

        'L^Tx=y から x を求める (後退代入)
        Dim x()  As Double = {0,0,0,0}
        backward_substitution(a,y,x)

        Console.WriteLine("X")
        disp_vector(x)
    End Sub

    '１次元配列の表示
    Private Sub disp_vector(ByVal row() As Double)
        For Each col As Double In row
            Console.Write(String.Format("{0,14:F10}{1}", col, vbTab))
        Next
        Console.WriteLine()
    End Sub

    '２次元配列の表示
    Private Sub disp_matrix(ByVal matrix(,) As Double)
        Dim i As Integer = 0
        For Each col As Double In matrix
            Console.Write(String.Format("{0,14:F10}{1}", col, vbTab))
            i += 1
            If i > N Then
                i = 0
                Console.WriteLine()
            End If
        Next
    End Sub

    '前進消去
    Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double)
        For pivot As Integer = 0 To N
            Dim s As Double = 0
            For col As Integer = 0 To (pivot - 1)
                s += a(pivot,col) * a(pivot,col)
            Next
            'ここで根号の中が負の値になると計算できない！
            a(pivot,pivot) = Math.Sqrt(a(pivot,pivot) - s)

            For row As Integer = (pivot + 1) To N
                s = 0
                For col As Integer = 0 To (pivot - 1)
                    s += a(row,col) * a(pivot,col)
                Next
                a(row,pivot) = (a(row,pivot) - s) / a(pivot,pivot)
                a(pivot,row) = a(row,pivot)
            Next
        Next
    End Sub

    '前進代入
    Private Sub forward_substitution(ByVal a(,) As Double, ByVal b() As Double, ByRef y() As Double)
        For row As Integer = 0 To N
            For col As Integer = 0 To row
                b(row) -= a(row,col) * y(col)
            Next
            y(row) = b(row) / a(row,row)
        Next
    End Sub

    '後退代入
    Private Sub backward_substitution(ByVal a(,) As Double, ByVal y() As Double, ByRef x() As Double)
        For row As Integer = N To 0 Step -1
            For col As Integer = N To (row + 1) Step - 1
                y(row) -= a(row,col) * x(col)
            Next
            x(row) = y(row) / a(row,row)
        Next
    End Sub
End Module

Z:\>vbc -nologo VB1006.vb

Z:\>VB1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

C#

using System;

public class CS1006
{
    private const int N = 4;

    // コレスキー法
    public static void Main()
    {
        double[,] a = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
        double[]  b = {34,68,96,125};

        Console.WriteLine("A");
        disp_matrix(a);
        Console.WriteLine("B");
        disp_vector(b);
        Console.WriteLine();

        // 前進消去
        forward_elimination(a,b);

        Console.WriteLine("LL^T");
        disp_matrix(a);

        // Ly=b から y を求める (前進代入)
        double[] y = {0,0,0,0};
        forward_substitution(a,b,y);

        Console.WriteLine("Y");
        disp_vector(y);

        // L^Tx=y から x を求める (後退代入)
        double[] x = {0,0,0,0};
        backward_substitution(a,y,x);

        Console.WriteLine("X");
        disp_vector(x);
    }
    // 前進消去
    private static void forward_elimination(double[,] a, double[] b)
    {
        for (int pivot = 0; pivot < N; pivot++)
        {
            double s = 0;
            for (int col = 0; col < pivot; col++)
                s += a[pivot,col] * a[pivot,col];
            // ここで根号の中が負の値になると計算できない！
            a[pivot,pivot] = Math.Sqrt(a[pivot,pivot] - s);

            for (int row = pivot + 1; row < N; row++)
            {
                s = 0;
                for (int col = 0; col < pivot; col++)
                    s += a[row,col] * a[pivot,col];
                a[row,pivot] =  (a[row,pivot] - s) / a[pivot,pivot];
                a[pivot,row] =  a[row,pivot];
            }
        }
    }
    // 前進代入
    private static void forward_substitution(double[,] a, double[] b, double[] y)
    {
        for (int row = 0; row < N; row++)
        {
            for (int col = 0; col < row; col++)
                b[row] -= a[row,col] * y[col];
            y[row] = b[row] / a[row,row];
        }
    }
    // 後退代入
    private static void backward_substitution(double[,] a, double[] y, double[] x)
    {
        for (int row = N - 1; row >= 0; row--)
        {
            for (int col = N - 1; col > row; col--)
                y[row] -= a[col,row] * x[col];
            x[row] = y[row] / a[row,row];
        }
    }
    // １次元配列を表示
    private static void disp_vector(double[] row)
    {
        foreach (double col in row)
            Console.Write(string.Format("{0,14:F10}\t", col));
        Console.WriteLine();
    }
    // ２次元配列を表示
    private static void disp_matrix(double[,] matrix)
    {
        int i = 0;
        foreach (double col in matrix)
        {
            Console.Write(string.Format("{0,14:F10}\t", col));
            if (++i >= N)
            {
                i = 0;
                Console.WriteLine();
            }
        }
    }
}

Z:\>csc -nologo CS1006.cs

Z:\>CS1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Java

public class Java1006 {

    private static final int N = 4;

    // コレスキー法
    public static void main(String []args) {
        double[][] a = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
        double[]   b = {34,68,96,125};

        System.out.println("A");
        disp_matrix(a);
        System.out.println("B");
        disp_vector(b);
        System.out.println();

        // 前進消去
        forward_elimination(a,b);

        System.out.println("LL^T");
        disp_matrix(a);

        // Ly=b から y を求める (前進代入)
        double y[] = {0,0,0,0};
        forward_substitution(a,b,y);

        System.out.println("Y");
        disp_vector(y);

        // L^Tx=y から x を求める (後退代入)
        double x[] = {0,0,0,0};
        backward_substitution(a,y,x);

        System.out.println("X");
        disp_vector(x);
    }
    // 前進消去
    private static void forward_elimination(double[][] a, double[] b) {
        for (int pivot = 0; pivot < N; pivot++) {
            double s = 0;
            for (int col = 0; col < pivot; col++)
                s += a[pivot][col] * a[pivot][col];
            // ここで根号の中が負の値になると計算できない！
            a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s);

            for (int row = pivot + 1; row < N; row++) {
                s = 0;
                for (int col = 0; col < pivot; col++)
                    s += a[row][col] * a[pivot][col];
                a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot];
                a[pivot][row] =  a[row][pivot];
            }
        }
    }
    // 前進代入
    private static void forward_substitution(double[][] a, double[] b, double[] y) {
        for (int row = 0; row < N; row++) {
            for (int col = 0; col < row; col++)
                b[row] -= a[row][col] * y[col];
            y[row] = b[row] / a[row][row];
        }
    }
    // 後退代入
    private static void backward_substitution(double[][] a, double[] y, double[] x) {
        for (int row = N - 1; row >= 0; row--) {
            for (int col = N - 1; col > row; col--)
                y[row] -= a[col][row] * x[col];
            x[row] = y[row] / a[row][row];
        }
    }
    // １次元配列を表示
    private static void disp_vector(double[] row) {
        for (double col: row)
            System.out.print(String.format("%14.10f\t", col));
        System.out.println();
    }
    // ２次元配列を表示
    private static void disp_matrix(double[][] matrix) {
        for (double[] row: matrix) {
            for (double col: row)
                System.out.print(String.format("%14.10f\t", col));
            System.out.println();
        }
    }
}

Z:\>javac Java1006.java

Z:\>java Java1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

C++

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const int N = 4;

// 前進消去
void forward_elimination(double a[N][N]);
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N]);
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// コレスキー法
int main()
{
    double a[N][N] = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
    double b[N]    = {34,68,96,125};

    cout << "A" << endl;
    disp_matrix(a);
    cout << "B" << endl;
    disp_vector(b);
    cout << endl;

    // 前進消去
    forward_elimination(a);

    cout << "LL^T" << endl;
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[N] = {0,0,0,0};
    forward_substitution(a, b, y);

    cout << "Y" << endl;
    disp_vector(y);

    // L^Tx=y から x を求める (後退代入)
    double x[N] = {0,0,0,0};
    backward_substitution(a, y, x);

    cout << "X" << endl;
    disp_vector(x);

    return 0;
}
// 前進消去
void forward_elimination(double a[N][N])
{
    for (int pivot = 0; pivot < N; pivot++)
    {
        double s = 0;
        for (int col = 0; col < pivot; col++)
            s += a[pivot][col] * a[pivot][col];
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = sqrt(a[pivot][pivot] - s);

        for (int row = pivot + 1; row < N; row++)
        {
            s = 0;
            for (int col = 0; col < pivot; col++)
                s += a[row][col] * a[pivot][col];
            a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot];
            a[pivot][row] =  a[row][pivot];
        }
    }
}
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N])
{
    for (int row = 0; row < N; row++)
    {
        for (int col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row] / a[row][row];
    }
}
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N])
{
    for (int row = N - 1; row >= 0; row--)
    {
        for (int col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col];
        x[row] = y[row] / a[row][row];
    }
}
// １次元配列を表示
void disp_vector(double row[N])
{
    for (int i = 0; i < N; i++)
        cout << setw(14) << fixed << setprecision(10) << row[i] << "\t";
    cout << endl;
}
// ２次元配列を表示
void disp_matrix(double matrix[N][N])
{
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
            cout << setw(14) << fixed << setprecision(10) << matrix[i][j] << "\t";
        cout << endl;
    }
}

Z:\>bcc32 -q CP1006.cpp
cp1006.cpp:

Z:\>CP1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Objective-C

#import <Foundation/Foundation.h>
#import <math.h>

const int N = 4;

// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N]);
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// コレスキー法
int main()
{
    double a[4][4] = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}};
    double b[4]    = {34,68,96,125};

    printf("A\n");
    disp_matrix(a);
    printf("B\n");
    disp_vector(b);
    printf("\n");

    // 前進消去
    forward_elimination(a,b);

    printf("LL^T\n");
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[4] = {0,0,0,0};
    forward_substitution(a,b,y);

    printf("Y\n");
    disp_vector(y);

    // L^Tx=y から x を求める (後退代入)
    double x[4] = {0,0,0,0};
    backward_substitution(a,y,x);

    printf("X\n");
    disp_vector(x);

    return 0;
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
    int pivot, row, col;
    for (pivot = 0; pivot < N; pivot++)
    {
        double s = 0;
        for (col = 0; col < pivot; col++)
            s += a[pivot][col] * a[pivot][col];
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = sqrt(a[pivot][pivot] - s);

        for (row = pivot + 1; row < N; row++)
        {
            s = 0;
            for (col = 0; col < pivot; col++)
                s += a[row][col] * a[pivot][col];
            a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot];
            a[pivot][row] =  a[row][pivot];
        }
    }
}
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N])
{
    int row, col;
    for (row = 0; row < N; row++)
    {
        for (col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row] / a[row][row];
    }
}
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N])
{
    int row, col;
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col];
        x[row] = y[row] / a[row][row];
    }
}
// １次元配列を表示
void disp_vector(double row[N])
{
    int i;
    for (i = 0; i < N; i++)
    {
        printf("%14.10f\t", row[i]);
    }
    printf("\n");
}
// ２次元配列を表示
void disp_matrix(double matrix[N][N])
{
    int i, j;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
            printf("%14.10f\t", matrix[i][j]);
        printf("\n");
    }
}

xxxxxx@yyyyyy /Z
$ gcc -o OC1006 OC1006.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

D

import std.stdio;
import std.math;

const int N = 4;

// コレスキー法
void main(string[] args)
{
    double[N][N] a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]];
    double[N]    b = [34,68,96,125];

    writefln("A");
    disp_matrix(a);
    writefln("B");
    disp_vector(b);
    writefln("");

    // 前進消去
    forward_elimination(a,b);

    writefln("LL^T");
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[N] = [0,0,0,0];
    forward_substitution(a,b,y);

    writefln("Y");
    disp_vector(y);

    // L^Tx=y から x を求める (後退代入)
    double x[N] = [0,0,0,0];
    backward_substitution(a,y,x);

    writefln("X");
    disp_vector(x);
}
// 前進消去
void forward_elimination(ref double[N][N] a, double[N] b)
{
    foreach (pivot; 0..N)
    {
        double s = 0;
        foreach (col; 0..pivot)
            s += a[pivot][col] * a[pivot][col];
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = sqrt(a[pivot][pivot] - s);

        foreach (row; (pivot + 1)..N)
        {
            s = a[row][pivot] / a[pivot][pivot];
            s = 0;
            foreach (col; 0..pivot)
                s += a[row][col] * a[pivot][col]; //
            a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot];
            a[pivot][row] =  a[row][pivot];
        }
    }
}
// 前進代入
void forward_substitution(double[N][N] a, double[N] b, ref double[N] y)
{
    foreach (row; 0..N)
    {
        foreach (col; 0..row)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row] / a[row][row];
    }
}
// 後退代入
void backward_substitution(double[N][N] a, double[N] y, ref double[N] x)
{
    foreach_reverse (row; 0..N)
    {
        foreach_reverse (col; (row+1)..N)
            y[row] -= a[col][row] * x[col];
        x[row] = y[row] / a[row][row];
    }
}
// １次元配列を表示
void disp_vector(double[] row)
{
    foreach(col; row)
        writef("%14.10f\t", col);
    writefln("");
}
// ２次元配列を表示
void disp_matrix(double[N][N] matrix)
{
    foreach(row; matrix)
    {
        foreach(col; row)
            writef("%14.10f\t", col);
        writefln("");
    }
}

Z:\>dmd D1006.d

Z:\>D1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Go

package main

import "fmt"
import "math"

const N = 4

// コレスキー法
func main() {
    var a [N][N]float64 = [N][N]float64{{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}}
    var b []float64     = []float64{34,68,96,125}

    fmt.Println("A")
    disp_matrix(a)
    fmt.Println("B")
    disp_vector(b)
    fmt.Println()

    // 前進消去
    forward_elimination(&a,b)

    fmt.Println("LL^T")
    disp_matrix(a)

    // Ly=b から y を求める (前進代入)
    var y []float64 = []float64{0,0,0,0}
    forward_substitution(a,b,y)

    fmt.Println("Y")
    disp_vector(y)

    // L^Tx=y から x を求める (後退代入)
    var x []float64 = []float64{0,0,0,0}
    backward_substitution(a,y,x)

    fmt.Println("X")
    disp_vector(x)
}
// 前進消去
func forward_elimination(a *[N][N]float64, b []float64) {
    for pivot := 0; pivot < N; pivot++ {
        var s = 0.0
        for col := 0; col < pivot; col++ {
            s += a[pivot][col] * a[pivot][col]
        }
        // ここで根号の中が負の値になると計算できない！
        a[pivot][pivot] = math.Sqrt(a[pivot][pivot] - s)

        for row := pivot + 1; row < N; row++ {
            s = a[row][pivot] / a[pivot][pivot]
            s = 0
            for col := 0; col < pivot; col++ {
                s += a[row][col] * a[pivot][col]
            }
            a[row][pivot] =  (a[row][pivot] - s) / a[pivot][pivot]
            a[pivot][row] =  a[row][pivot]
        }
    }
}
// 前進代入
func forward_substitution(a [N][N]float64, b []float64, y []float64) {
    for row := 0; row < N; row++ {
        for col := 0; col < row; col++ {
            b[row] -= a[row][col] * y[col]
        }
        y[row] = b[row] / a[row][row]
    }
}
// 後退代入
func backward_substitution(a [N][N]float64, y []float64, x []float64) {
    for row := N - 1; row >= 0; row-- {
        for col := N - 1; col > row; col-- {
            y[row] -= a[col][row] * x[col]
        }
        x[row] = y[row] / a[row][row]
    }
}
// １次元配列を表示
func disp_vector(row[]float64) {
    for _, col := range row {
        fmt.Printf("%14.10f\t", col)
    }
    fmt.Println("")
}
// ２次元配列を表示
func disp_matrix(matrix[N][N]float64) {
    for _, row := range matrix {
        for _, col := range row {
            fmt.Printf("%14.10f\t", col)
        }
        fmt.Println("")
    }
}

Z:\>8g GO1006.go

Z:\>8l -o GO1006.exe GO1006.8

Z:\>GO1006
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Scala

object Scala1006 {
    val N = 3

    def main(args: Array[String]) {
        var a:Array[Array[Double]] = Array(Array(5,2,3,4), Array(2,10,6,7), Array(3,6,15,9), Array(4,7,9,20))
        var b:Array[Double] = Array(34,68,96,125)

        println("A")
        disp_matrix(a)
        println("B")
        disp_vector(b)
        println()

        // 前進消去
        forward_elimination(a,b)

        println("LL^T")
        disp_matrix(a)

        // Ly=b から y を求める (前進代入)
        var y:Array[Double] = Array(0,0,0,0)
        forward_substitution(a,b,y)

        println("Y")
        disp_vector(y)

        // L^Tx=y から x を求める (後退代入)
        var x:Array[Double] = Array(0,0,0,0)
        backward_substitution(a,y,x)

        println("X")
        disp_vector(x)
    }
    // 前進消去
    def forward_elimination(a:Array[Array[Double]], b:Array[Double]) {
        for (pivot <- 0 to N) {
            var s:Double = 0;
            for (col <- 0 to pivot - 1)
                s += a(pivot)(col) * a(pivot)(col)
            // ここで根号の中が負の値になると計算できない！
            a(pivot)(pivot) = Math.sqrt(a(pivot)(pivot) - s)

            for (row <- pivot + 1 to N) {
                s = 0
                for (col <- 0 to pivot - 1)
                    s += a(row)(col) * a(pivot)(col)
                a(row)(pivot) =  (a(row)(pivot) - s) / a(pivot)(pivot)
                a(pivot)(row) =  a(row)(pivot)
            }
        }
    }
    // 前進代入
    def forward_substitution(a:Array[Array[Double]], b:Array[Double], y:Array[Double]) {
        for (row <- 0 to N) {
            for (col <- 0 to row - 1)
                b(row) -= a(row)(col) * y(col)
            y(row) = b(row) / a(row)(row)
        }
    }
    // 後退代入
    def backward_substitution(a:Array[Array[Double]], y:Array[Double], x:Array[Double]) {
        for (row <- (0 to N).reverse) {
            for (col <- (row + 1 to N).reverse)
                y(row) -= a(row)(col) * x(col)
            x(row) = y(row) / a(row)(row)
        }
    }

    // １次元配列を表示
    def disp_vector(row:Array[Double]) = {
        row.foreach { col =>
            print("%14.10f\t".format(col))
        }
        println()
    }
    // ２次元配列を表示
    def disp_matrix(matrix:Array[Array[Double]]) = {
        matrix.foreach { row =>
            row.foreach { col =>
                print("%14.10f\t".format(col))
            }
            println()
        }
    }
}

Z:\>scala Scala1006.scala
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

F#

module Fs1006

open System

let N = 3

// １次元配列を表示
let disp_vector (row:float[]) =
    row
    |> Array.iter (fun x -> printf "%14.10f" x)
    printfn ""

// ２次元配列を表示
let disp_matrix (matrix:float[][]) =
    matrix
    |> Array.iter
        (fun row ->
            row
            |> Array.iter (fun col -> printf "%14.10f" col)
            printfn ""
        )

// 前進消去
let forward_elimination (a:float[][]) (b:float[]) =
    for pivot in [0..N] do
        let mutable s:float = 0.0
        for col in [0..(pivot - 1)] do
            s <- s + a.[pivot].[col] * a.[pivot].[col]
        // ここで根号の中が負の値になると計算できない！
        a.[pivot].[pivot] <- Math.Sqrt(a.[pivot].[pivot] - s)


        for row in [pivot + 1 .. N] do
            s <- 0.0
            for col in [0..(pivot-1)] do
                s <- s + a.[row].[col] * a.[pivot].[col]
            a.[row].[pivot] <- (a.[row].[pivot] - s) / a.[pivot].[pivot]
            a.[pivot].[row] <- a.[row].[pivot]

// 前進代入
let forward_substitution (a:float[][]) (b:float[]) (y:float[]) =
    for row in [0..N] do
        for col in 0..(row - 1) do
            b.[row] <- b.[row] - a.[row].[col] * y.[col]
        y.[row] <- b.[row] / a.[row].[row]

// 後退代入
let backward_substitution (a:float[][]) (y:float[]) (x:float[]) =
    for row in (List.rev [0..N]) do
        for col in (List.rev [row + 1 .. N]) do
            y.[row] <- y.[row] - a.[row].[col] * x.[col]
        x.[row] <- y.[row] / a.[row].[row]

let mutable a:float[][] = [| [|5.0;2.0;3.0;4.0 |]; [|2.0;10.0;6.0;7.0 |]; [| 3.0;6.0;15.0;9.0 |]; [|4.0;7.0;9.0;20.0 |] |]
let mutable b:float[] = [| 34.0;68.0;96.0;125.0 |]

printfn "A"
disp_matrix a
printfn "B"
disp_vector b
printfn ""

// 前進消去
forward_elimination a b

printfn "LL^T"
disp_matrix a

// Ly=b から y を求める (前進代入)
let mutable y:float[] = [| 0.0;0.0;0.0;0.0 |]
forward_substitution a b y

printfn "Y"
disp_vector y

// L^Tx=y から x を求める (後退代入)
let mutable x:float[] = [| 0.0;0.0;0.0;0.0 |]
backward_substitution a y x

printfn "X"
disp_vector x

exit 0

Z:\>fsi  --nologo --quiet Fs1006.fs
A
  5.0000000000  2.0000000000  3.0000000000  4.0000000000
  2.0000000000 10.0000000000  6.0000000000  7.0000000000
  3.0000000000  6.0000000000 15.0000000000  9.0000000000
  4.0000000000  7.0000000000  9.0000000000 20.0000000000
B
 34.0000000000 68.0000000000 96.0000000000125.0000000000

LL^T
  2.2360679775  0.8944271910  1.3416407865  1.7888543820
  0.8944271910  3.0331501776  1.5825131361  1.7803272782
  1.3416407865  1.5825131361  3.2704207946  1.1566122322
  1.7888543820  1.7803272782  1.1566122322  3.5060922587
Y
 15.2052622470 17.9351488764 14.4377113129 14.0243690350
X
  1.0000000000  2.0000000000  3.0000000000  4.0000000000

Clojure

(def N 4)

(def a [[5.0 2.0 3.0 4.0] [2.0 10.0 6.0 7.0] [3.0 6.0 15.0 9.0] [4.0 7.0 9.0 20.0]])
(def b  [34.0 68.0 96.0 125.0])

;１次元配列を表示
(defn disp_vector [row]
    (doseq [col row]
        (print (format "%14.10f\t" col)))
    (println))

;２次元配列を表示
(defn disp_matrix [matrix]
    (doseq [row matrix]
        (doseq [col row]
            (print (format "%14.10f\t" col)))
        (println)))

;前進消去
(defn forward_elim_loop [pivot row a a_sqrt]
    (def a_map (map
        (fn [a_row a_piv] (* a_row a_piv))
            (take pivot (nth a row)) (take pivot (nth a pivot))))
    (def s (apply + a_map))
    (def a_col (/ (- (nth (nth a row) pivot) s) a_sqrt))

    (def a_row (concat (take pivot (nth a row  )) (cons a_col (drop (inc pivot) (nth a row)))))
    (def a_piv (concat (take row   (nth a pivot)) (cons a_col (drop (inc row)   (nth a pivot)))))

    (def a1 (concat (take pivot a ) (cons (vec a_piv) (drop (inc pivot) a))))
    (def a2 (concat (take row   a1) (cons (vec a_row) (drop (inc row)   a))))

    (if (< row (dec N))
        (forward_elim_loop pivot (inc row) a2 a_sqrt)
        a2))

(defn forward_elimination [pivot a]
    (def a_map (map (fn [col] (* col col)) (take pivot (nth a pivot))))
    (def s (apply + a_map))

    ;ここで根号の中が負の値になると計算できない！
    (def a_sqrt (. Math sqrt (- (nth (nth a pivot) pivot) s)))
    (def a2 (forward_elim_loop pivot pivot a a_sqrt))

    (if (< pivot (dec N))
        (forward_elimination (inc pivot) a2)
        a2))

;前進代入
(defn forward_substitution [row a b]
    (def a_map (map
        (fn [a_col b_col] (* a_col b_col))
            (take row (nth a row)) b))
    (def s (apply + a_map))
    (def y (/ (- (nth b row) s) (nth (nth a row) row)))
    (def b2 (concat (take row b) (cons y (drop (inc row) b))))
    (if (< row (dec N))
        (cons y (forward_substitution (inc row) a b2))
        (cons y [])))

;後退代入
(defn backward_substitution [row a b]
    (def a_map (map
        (fn [a_col b_col] (* a_col b_col))
            (drop (inc row) (nth a row)) (drop (inc row) b)))
    (def s (apply + a_map))
    (def x (/ (- (nth b row) s) (nth (nth a row) row)))
    (def b2 (concat (take row b) (cons x (drop (inc row) b))))
    (if (> row 0)
        (cons x (backward_substitution (dec row) a b2))
        (cons x [])))

(println "A")
(disp_matrix a)
(println "B")
(disp_vector b)
(println "")

;前進消去
(def a2 (forward_elimination 0 a))

(println "LL^T")
(disp_matrix a2)

;Ly=b から y を求める (前進代入)
(def y (forward_substitution 0 a2 b))

(println "Y")
(disp_vector y)

;L^Tx=y から x を求める (後退代入)
(def x (backward_substitution (dec N) a2 y))

(println "X")
(disp_vector (reverse x))

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj1006.clj
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

Haskell

import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- ２次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- １次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->Double->[[Double]]
forward_elim_loop pivot row a a_sqrt =
    let
        a_zip = zip (take pivot (a!!row)) (take pivot (a!!pivot))
        s     = sum $ map(\(a_row, a_pivot) -> a_row * a_pivot) a_zip
        a_col = (a!!row!!pivot - s) / a_sqrt

        a_row = (take pivot (a!!row)  ) ++ (a_col:(drop (pivot+1) (a!!row)))
        a_piv = (take row   (a!!pivot)) ++ (a_col:(drop (row+1)   (a!!pivot)))

        a1 = (take pivot a ) ++ (a_piv:(drop (pivot+1) a))
        a2 = (take row   a1) ++ (a_row:(drop (row+1)   a))
    in
        if row < (n - 1)
            then forward_elim_loop pivot (row+1) a2 a_sqrt
            else a2

forward_elimination::Int->[[Double]]->[[Double]]
forward_elimination pivot a =
    let
        s      = sum $ map(\col -> col * col) $ take pivot (a!!pivot)
        -- ここで根号の中が負の値になると計算できない！
        a_sqrt = (sqrt (a!!pivot!!pivot - s))
        a2     = forward_elim_loop pivot pivot a a_sqrt
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2
            else
                a2

-- 前進代入
forward_substitution::Int->[[Double]]->[Double]->[Double]
forward_substitution row a b =
    let
        a_zip = zip (take row (a!!row)) b
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        y  = (b!!row - s) / a!!row!!row
        b2 = (take row b) ++ (y:(drop (row+1) b))
    in
        if row < (n - 1)
            then y:(forward_substitution (row+1) a b2)
            else y:[]

-- 後退代入
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20::Double]]
    let b  = [34,68,96,125::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    -- 前進消去
    let a2 = forward_elimination 0 a

    putStrLn "LL^T"
    disp_matrix a2

    -- Ly=b から y を求める (前進代入)
    let y = forward_substitution 0 a2 b

    putStrLn "Y"
    disp_vector y

    -- L^Tx=y から x を求める (後退代入)
    let x = backward_substitution (n-1) a2 y

    putStrLn "X"
    disp_vector (reverse x)

Z:\>runghc Hs1006.hs
A
  5.0000000000    2.0000000000    3.0000000000    4.0000000000
  2.0000000000   10.0000000000    6.0000000000    7.0000000000
  3.0000000000    6.0000000000   15.0000000000    9.0000000000
  4.0000000000    7.0000000000    9.0000000000   20.0000000000
B
 34.0000000000   68.0000000000   96.0000000000  125.0000000000

LL^T
  2.2360679775    0.8944271910    1.3416407865    1.7888543820
  0.8944271910    3.0331501776    1.5825131361    1.7803272782
  1.3416407865    1.5825131361    3.2704207946    1.1566122322
  1.7888543820    1.7803272782    1.1566122322    3.5060922587
Y
 15.2052622470   17.9351488764   14.4377113129   14.0243690350
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

VBScript
JScript
PowerShell
Perl
PHP
Python
Ruby
Groovy

Pascal
Ada
VB.NET

C#
Java
C++
Objective-C
D
Go

Scala
F#
Clojure
Haskell